
Calhoun 

iniQiuiic^iul Ar{hiv« of tilt Mil vdl Poii^roduiit School 


Calhoun: The NPS Institutional Archive 
□Space Repository 



Theses and Dissertations 


1. Thesis and Dissertation Collection, all items 


2014-12 

Galerkin optimal control 

Boucher, Randy 

Monterey, California: Naval Postgraduate School 


http://hdl.handle.net/10945/44526 


This publication is a work of the U.S. Government as defined in Title 17, United 
States Code, Section 101. Copyright protection is not available for this work in the 
United States. 

Downloaded from NPS Archive: Calhoun 



DUDLEY 

KNOX 

LIBRARY 


htt p://w ww. n ps. e du/l ib ra ry 


Caflwuo is the Naval Postgraduate School's public access digital repository for 
research mate rials and institutiional publicatkins created by the NPS community. 
Calhoun is named for Professor of Mathematics Guy K. Caftiouo, NPS's first 
appointed — and published — schoteily author. 

Dudley Knox Library / Naval Postgraduate School 
411 Dyer Road / 1 Univefsity Circle 
Monterey, California USA 93943 







NAVAL 

POSTGRADUATE 

SCHOOL 

MONTEREY, CALIFORNIA 

DISSERTATION 


GALERKIN OPTIMAL CONTROL 

by 

Randy Boucher 
December 2014 


Dissertation Supervisor: 


Wei Kang 


Approved for public release; distribution is unlimited 




THIS PAGE INTENTIONALLY LEET BLANK 



REPORT DOCUMENTATION PAGE 

Form Approved 0MB No. 0704-0188 

Public reporting burden for this collection of information is estimated to average 1 hour per response, including the time for reviewing instruction, searching existing data sources, gathering and 
maintaining the data needed, and completing and reviewing the collection of information. Send comments regarding this burden estimate or any other aspect of this collection of information, 
including suggestions for reducing this burden, to Washington headquarters Services, Directorate for Information Operations and Reports, 1215 Jefferson Davis Highway, Suite 1204, Arlington, 

VA 22202-4302, and to the Office of Management and Budget, Paperwork Reduction Project (0704-0188) Washington DC 20503. 

1 . AGENCY USE ONLY (Leave blank) 2. REPORT DATE 3. REPORT TYPE AND DATES COVERED 

December 2014 Dissertation - October 13 - December 14 

4. TITLE AND SUBTITLE: GALERKIN OPTIMAL CONTROL 

5. FUNDING NUMBERS 

6. AUTHOR(S): Randy Boucher 

7. PERFORMING ORGANIZATION NAME(S) AND ADDRESS(ES) 

Naval Postgraduate School 

Monterey, CA 93943-5000 

8. PERFORMING ORGANIZATION 
REPORT NUMBER 

9.SPONSORINGMONITORING AGENCY NAME(S) AND ADDRESS(ES) 

10. SPONSORING/MONITORING 
AGENCY REPORT NUMBER 

11. SUPPLEMENTARY NOTES: The views expressed in this thesis are those of the author and do not reflect the 
official policy or position of the Department of Defense or the U.S. Government. 

12a. DISTRIBUTION / AVAILABILITY STATEMENT 

Approved for public release; distribution is unlimited 

12b. DISTRIBUTION CODE 

A 

13. ABSTRACT (maximum 200 words) 

A Galerkin-based family of numerical formulations is presented for solving nonlinear optimal control problems. This 
dissertation introduces a family of direct methods that calculate optimal trajectories by discretizing the system dy¬ 
namics using Galerkin numerical techniques and approximate the cost function with Gaussian quadrature. In this 
numerical approach, the analysis is based on L^-norms. An important result in the theoretical foundation is that the 
feasibility and consistency theorems are proved for problems with continuous and/or piecewise continuous controls. 
Galerkin methods may be formulated in a number of ways that allow for efficiency and/or improved accuracy while 
solving a wide range of optimal control problems with a variety of state and control constraints. Numerical formula¬ 
tions using Lagrangian and Legendre test functions are derived. One formulation allows for a weak enforcement of 
boundary conditions, which imposes end conditions only up to the accuracy of the numerical approximation itself. 
Additionally, the multi-scale formulation can reduce the dimension of multi-scale optimal control problems, those in 
which the states and controls evolve on different timescales. Finally, numerical examples are shown to demonstrate 
the versatile nature of Galerkin optimal control. 

14. SUBJECT TERMS 

Galerkin, Pseudospectral, Optimal Control, Constrainted Optimization 

15. NUMBER 

OF PAGES 

217 

16. PRICE CODE 

17. SECURITY CLASSIFICATION 
OF REPORT 

Unclassified 

18. SECURITY CLASSIFICATION 

OF THIS PAGE 

Unclassified 

19. SECURITYCLASSIFICATION 

OF ABSTRACT 

Unclassified 

20. LIMITATION OF 
ABSTRACT 

uu 


NSN 7540-01-280-5500 


Standard Form 298 (Rev. 2-89) 
Prescribed by ANSI Std. 239-18 


1 




























THIS PAGE INTENTIONALLY LEET BLANK 


11 



Approved for public release; distribution is unlimited 
GALERKIN OPTIMAL CONTROL 


Randy Boucher 

Lieutenant Colonel, United States Army 
B.S., Boston University, 1997 
M.S., University of Washington, 2006 

Submitted in partial fulfillment of the 
requirements for the degree of 

DOCTOR OF PHILOSOPHY IN 
APPLIED MATHEMATICS 

from the 

NAVAL POSTGRADUATE SCHOOL 
December 2014 


Author: 

Randy Boucher 


Approved By: 

Wei Kang 

Professor 

Department of Appl. Math. 
Dissertation Supervisor 

Chris L. Frenzen 

Professor 

Department of Appl. Math. 


Francis X. Giraldo 

Professor 

Department of Appl. Math. 

Arthur J. Kroner 

Professor 

Department of Appl. Math. 


I. Michael Ross 

Professor 

Department of Mech. & Aero. 


Approved By: 

Craig Rasmussen, Professor and Chair, Department of Appl. Math. 

Approved By: 

Doug Moses, Vice Provost for Academic Affairs 



THIS PAGE INTENTIONALLY LEET BLANK 


IV 



ABSTRACT 


A Galerkin-based family of numerical formulations is presented for solving nonlinear op¬ 
timal control problems. This dissertation introduces a family of direct methods that cal¬ 
culate optimal trajectories by discretizing the system dynamics using Galerkin numerical 
techniques and approximate the cost function with Gaussian quadrature. In this numer¬ 
ical approach, the analysis is based on L^-norms. An important result in the theoretical 
foundation is that the feasibility and consistency theorems are proved for problems with 
continuous and/or piecewise continuous controls. Galerkin methods may be formulated in 
a number of ways that allow for efficiency and/or improved accuracy while solving a wide 
range of optimal control problems with a variety of state and control constraints. Numerical 
formulations using Lagrangian and Legendre test functions are derived. One formulation 
allows for a weak enforcement of boundary conditions, which imposes end conditions only 
up to the accuracy of the numerical approximation itself. Additionally, the multi-scale 
formulation can reduce the dimension of multi-scale optimal control problems, those in 
which the states and controls evolve on different timescales. Finally, numerical examples 
are shown to demonstrate the versatile nature of Galerkin optimal control. 


V 



THIS PAGE INTENTIONALLY LEET BLANK 


VI 



TABLE OF CONTENTS 


1 INTRODUCTION . 1 

2 MATHEMATICAL BACKGROUND. 3 

2.1 Optimal Control. 3 

2.1.1 The Optimal Control Problem. 4 

2.1.2 Pontryagin’s Minimum Principle. 5 

2.1.2.1 Calculus of Variations. 5 

2.1.2.2 Necessary Conditions . 6 

2.1.2.3 Pontryagin’s Minimum Principle. 7 

2.1.3 Numerical Methods for Optimal Control . 9 

2.2 Interpolation and Numerical Integration. 11 

2.2.1 Interpolation. 11 

2.2.1.1 Modal Interpolation. 11 

2.2.1.2 Nodal Interpolation . 14 

2.2.1.3 Transformations between grids. 16 

2.2.1.4 Interpolation Quality. 18 

2.2.2 Numerical Integration. 25 

2.2.2.1 Gaussian Quadrature. 27 

2 . 22.2 Gaussian Quadrature Accuracy. 28 

2.3 Numerical Solutions to Differential Equations. 31 

2.3.1 Spectral Methods . 32 

2.3.1.1 Collocation. 32 

2.3.1.2 Galerkin Numerical Methods. 35 

2.3.2 Spectral Element Methods . 43 

2.3.2.1 Continuous Galerkin. 43 

2.3.2.2 Discontinuous Galerkin. 48 

3 MOTIVATION FOR GALERKIN OPTIMAL CONTROL. 55 

3.1 Problem B (Bolza Problem). 55 

3.2 Legendre Pseudospectral Method. 56 

3.3 Problem B (Multi-scale Bolza Problem). 62 

3.4 A Modified Legendre PS Method for Multi-scale Problems.63 

3.5 Jackson’s Theorem. 67 

3.6 Motivation for the Weak Integral Formulation. 70 

vii 



































4 GENERAL GALERKIN OPTIMAL CONTROL FORMULATION.77 

4.1 Method for Approximation. 77 

4.2 Computation Strategy. 82 

4.2.1 Computation Strategy for GOCM-S. 82 

4.2.2 Computation Strategy for GOCM-S. 83 

4.3 Feasibility of Solutions . 83 

4.3.1 Feasibility of GOCM-S . 89 

4.3.2 Feasibility of GOCM-S . 92 

4.4 Consistency of Solutions .. 96 

4.4.1 Consistency of GOCM-S . 96 

4.4.2 Consistency of GOCM-S .103 

5 ALTERNATIVE FORMS OF GALERKIN OPTIMAL CONTROL.109 

5.1 Galerkin optimal control with Weak Boundary Condition Enforcement 109 

5.1.1 Method for Approximation.110 

5.1.2 Computation Strategy.114 

5.1.2.1 Computation Strategy for GOCM-W.114 

5.1.2.2 Computation Strategy for GOCM-W.115 

5.1.3 Feasibility of Solutions.116 

5.1.3.1 Feasibility of GOCM-W.116 

5.1.3.2 Feasibility of GOCM-W.119 

5.2 Galerkin Optimal Control with Continuous Element-based Galerkin . 122 

5.2.1 Method for Approximation.122 

5.2.2 Computation Strategy.126 

5.3 GOCM with Discontinuous Element-based Galerkin.127 

5.3.1 Method for Approximation.127 

5.3.2 Computation Strategy.131 

5.4 Galerkin Optimal Control for Multi-scale Problems.132 

5.4.1 Method for Approximation.132 

5.4.2 Computation Strategy.135 

5.5 Modifications to Galerkin Optimal Control.136 

5.5.1 Over-Integration of the RHS Vector.136 

5.5.2 Galerkin Optimal Control with LG and LGR/F-LGR Quadra¬ 
ture Points.138 

5.5.2.1 Method for Approximation for Galerkin Optimal Con¬ 
trol with LGR/F-LGR Nodes.139 

5.5.2.2 Method for Approximation for Galerkin Optimal Con¬ 
trol with LG Nodes.142 


viii 



































6 GALERKIN OPTIMAL CONTROL WITH LEGENDRE POLYNOMIAL 

TEST FUNCTIONS .147 

6.1 Methods for Approximation.147 

6.2 Computation Strategy.150 

6.2.1 Computation Strategy for GOCM-L.150 

6.2.2 Computation Strategy for GOCM-L.151 

6.3 Feasibility of Solutions .151 

6.4 Consistency of Solutions .154 

7 EXAMPLE PROBLEMS.159 

7.1 Example 7.1: Nonlinear Two-Dimensional Problem with Fixed Initial 

Conditions.159 

7.2 Example 7.2: Nonlinear Two-Dimensional Problem with Fixed Bound¬ 
ary Conditions.164 

7.3 Example 7.3: Two-Dimensional Bang-Bang Control Problem with Fixed 

Boundary Conditions.168 

7.4 Example 7.4: Two-Dimensional Problem with Fixed Boundary Condi¬ 
tions .172 

7.5 Example 7.5: Two-Dimensional Multi-scale Problem with Fixed Bound¬ 
ary Conditions.176 

7.6 Example 7.6: Nonlinear Two-Dimensional Multi-scale Problem with 

Fixed Initial Conditions.179 

7.7 Summary.183 

8 CONCLUSIONS AND AREAS FOR FUTURE RESEARCH.185 

8.1 Dissertation Summary .185 

8.2 Future Work.186 

APPENDIX A. NORMS AND FUNCTIONAL SPACES.187 

APPENDIX B. DISCRETE NORMS.189 

APPENDIX C. SOBOLEV SPACES.191 

LIST OF REFERENCES.193 

INITIAL DISTRIBUTION LIST.201 


IX 


























THIS PAGE INTENTIONALLY LEET BLANK 


X 



LIST OF FIGURES 


Figure 1. 
Figure 2. 
Figure 3. 
Figure 4. 
Figure 5. 
Figure 6. 
Figure 7. 
Figure 8. 
Figure 9. 
Figure 10. 
Figure 11. 
Figure 12. 
Figure 13. 
Figure 14. 

Figure 15. 
Figure 16. 
Figure 17. 
Figure 18. 
Figure 19. 
Figure 20. 
Figure 21. 
Figure 22. 
Figure 23. 
Figure 24. 
Figure 25. 
Figure 26. 
Figure 27. 
Figure 28. 
Figure 29. 
Figure 30. 
Figure 31. 


Legendre polynomials. 12 

Lagrange polynomials defined on equi-spaced grid. 15 

LG points. 19 

LGL points. 20 

Lagrange polynomials defined on LGL grid. 20 

LGR points. 21 

Plot of f{t) = cos(37rt). 22 

Interpolation of f{t) = cos(37rt) with 10-th order LGL points.23 

Interpolation of f{t) = cos(37rt) with 30-th order LGL points.24 

Comparison of interpolation errors for various orders of iV.25 

Comparison of quadrature errors for various orders of iV.30 

Exact solution and Legendre PS method approximation for Example 3.1. 60 

Eegendre spectral coefficients for Example 3.1.61 

Exact solution and multi-scale Eegendre PS method approximation for 

Example 3.1. 66 

Exact solution and GOCM-MS approximation for Example 3.1. 74 

Exact solution and GOCM-W approximation for Example 7.1.161 


State Xi approximation error vs. polynomial order, N, for Example 7.1. 162 
State X 2 approximation error vs. polynomial order, N, for Example 7.1. 163 
Control u approximation error vs. polynomial order, N, for Example 7.1. 163 

Exact solution and GOCM-W approximation for Example 7.2.165 

State Xi approximation error vs. polynomial order, N, for Example 7.2. 166 
State X 2 approximation error vs. polynomial order, N, for Example 7.2. 167 
Control u approximation error vs. polynomial order, N, for Example 7.2. 167 


Exact solution and GOCM-S approximation for Example 7.3.169 

Exact solution and GOCM-W approximation for Example 7.3.170 

Exact solution and GOCM-CG approximation for Example 7.3.171 

Exact solution and GOCM-DG approximation for Example 7.3.172 

Exact solution and GOCM-W approximation for Example 7.4.173 


State Xi approximation error vs. polynomial order, N, for Example 7.4. 174 
State X 2 approximation error vs. polynomial order, N, for Example 7.4. 175 
Control u approximation error vs. polynomial order, N, for Example 7.4. 175 


XI 
























Figure 32. Exact solution and GOCM-MS approximation with = 10, = 

50 and Nu = 5 for Example 7.5.177 

Eigure 33. Exact solution and GOCM-MS approximation with = 3, = 43 

and Nu = 1 for Example 7.5.178 

Eigure 34. Exact state solutions and GOCM-W approximation for Example 7.6. . . 180 
Eigure 35. Exact control solution and GOCM-W approximation for Example 7.6. . 180 

Eigure 36. Eegendre spectral coefficients for Example 7.6.181 

Eigure 37. Exact state solutions and GOCM-MS approximation for Example 7.6. . 182 
Eigure 38. Exact control solution and GOCM-MS approximation for Example 7.6. . 183 






ACKNOWLEDGEMENTS 


I would like to begin by thanking my wife, Jill. You are my compass and my 
strength. Thank you for keeping me anchored to all things that truly matter in life. Thanks 
also to my beautiful children, Leah, Dominic and Stella, who greet me with love and a 
smile at the end of each day. 

I would also like to acknowledge my wonderful team of researchers. Thanks espe¬ 
cially to Dr. Wei Kang for your patience and direction. You always made me feel like a 
colleague although you are one of many giants on whose shoulders I stand. Thanks also to 
the members of my committee: Dr. Frank Giraldo, Dr. Chris Frenzen, Dr. Mike Ross and 
Dr. Arthur Krener. Each of you planted great seeds of knowledge within me for which I 
am forever grateful. 

I have encountered many great people along this journey. Thank you to all of my 
friends at NFS for lending an ear and raising a pint! 



THIS PAGE INTENTIONALLY LEET BLANK 


XIV 



CHAPTER 1: 
INTRODUCTION 


The last two decades have proven to be a time of active research for numerical meth¬ 
ods for optimal control. Particularly, direct collocation methods, such as pseudospectral 
(PS) methods, have received much attention [1-8]. PS methods produce accurate solutions 
on a wide variety of optimal control problems. Two recent highlights are the success¬ 
ful use of the Legendre PS method for the first ever zero-propellant attitude maneuver of 
the International Space Station [4] and the first ever minimum-time rotational maneuver 
performed in orbit by a NASA space telescope called TRACE [8]. In the Legendre PS 
method [I, 2, 5, 6, 9, 10], the problem is discretized at the Legendre-Gauss-Lobatto (LGL) 
points, Legendre-Gauss-Radau (LGR) points or Legendre-Gauss (LG) points. The states 
are approximated with globally interpolating Lagrange polynomials and the cost function 
is typically approximated using Gaussian quadrature rule. Other variants of the PS method 
include the Chebyshev PS method [11], the PS knotting method [12] and the Bellman 
method [13]. The Legendre PS method will be outlined in Chapter 3 of this dissertation, 
preceded by a review of mathematical topics in Chapter 2. 

While PS methods have shown to be good all-round methods for solving nonlinear 
optimal control problems, approximating the derivative of a function using a standard PS 
differentiation matrix (such as the Legendre PS differentiation matrix) may introduce errors 
into the approximation. Chapter 3 highlights this issue with the use of Jackson’s Theorem. 
Additionally, Chapter 3 will motivate the use of the weak integral formulation to approx¬ 
imate the system’s dynamics. This leads to the creation of a family of Galerkin-based 
formulations called, “Galerkin optimal control.” 

The family of methods proposed in this dissertation are derived from Galerkin nu¬ 
merical techniques that have been developed for numerical solutions to differential equa¬ 
tions since the early 1970s [14—16]. In addition to the family of Galerkin optimal control 
formulations that are presented, this dissertation highlights important theorems that prove 
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method feasibility and consistency for problems with continuous and/or piecewise contin¬ 
uous controls. 

The base Galerkin optimal control method is outlined in Chapter 4, where fea¬ 
sibility and convergence theorems are presented. Chapter 5 presents a review of addi¬ 
tional Galerkin-based formulations and strategies such as the use element-based Galerkin 
techniques and a multi-scale approach. Lastly, modifications to the method such as over¬ 
integration and the use of various quadrature rules are offered to improve computational 
efficiency and/or increase accuracy of the solutions. 

The remainder of the dissertation is organized as follows: Chapter 6 presents a 
Petrov-Galerkin optimal control approach to discretizing the system dynamics; in place of 
Lagrange polynomial test functions integrated into the base formulation, a set of Legen¬ 
dre polynomials are used. Improved feasibility and convergence theorems are presented. 
Chapter 7 demonstrates the versatile nature of the Galerkin optimal control formulations by 
considering a number of example problems. Lastly, Chapter 8, highlights the potential for 
Galerkin optimal control in solving a wide range of real-world optimal control problems 
with a variety of state and control constraints. Additionally, Chapter 8 discusses areas of 
future research. 
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CHAPTER 2: 

MATHEMATICAL BACKGROUND 


2.1. Optimal Control 

Optimal control has a rich history that dates back to 1696, when Johann Bernoulli 
posed the bachristochrone problem in the Acta Eruditorum to [17, 18] ""the most astute 
mathematicians of the world." The bachristochrone problem was the following: 

If in the vertical plane two points A and B are given, then it is required to 
specify the orbit AMB of the movable point M, along which it, starting from A, 
and under the influence of its own weight, arrives at B in the shortest possible 
time. [19] 

In addition to Johann Bernoulli, other mathematical giants living in Europe at this time, 
such as Newton, Eeibniz and Johann’s brother, Jacob Bernoulli [19] (all considered “Men 
of Mathematics” by Bell [20]), solved the bachristochrone problem. Eater, Euler invented a 
method for solving such problems (with mathematical underpinnings created by Eagrange), 
known today as the foundations of the calculus of variations. The standard calculus of 
variations problems is of the form [21] 

minimize J[r/(-)] = / F{t,y{t),y{t))dt, (2.1) 

Jto 

subject to yfo) = r/o and y{tf) = yf, (2.2) 

where J acts on a set of functions and is called di functional. Notice that problem (2. l)-(2.2) 
may be written in the equivalent optimal control problem form 

minimize J[a;(-),«(■)] = / F{x{t),u{t))dt, (2.3) 

7 to 

subject to a;(to) = [to,yof', x{tf) = [tf,yff and x(t) = te [to,tf], (2.4) 
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by renaming the variables t and y as t = yi and y = y2, and creating the new vector 
X = [yuy2r. 

Over 250 years later, after many periods of active research in the field of calculus 
of variations, the Russian mathematician Lev Semenovich Pontryagin made a giant leap 
forward. In 1956, Pontryagin and his group established the optimal control theory [22, 
23]. In contrast to standard calculus of variations problems of the form (2.1)-(2.2), or 
equivalent form (2.3)-(2.4), it was shown that optimal control theory was well suited to 
handle discontinuous solutions, u{t). Additionally, Pontryagin established that problems 
of optimal control involved the minimization of a functional over a set of function pairs, 
t H-)■ {x, u) G X subject to the dynamical constraint 

x{t) = f{x{t),u{t)), 

where / : x —)■ and u{t) is a control function. It was soon realized that this 

new theory of optimal control was well suited to solve many complex problems (that the 
calculus of variations could not). Over the last half century, optimal control theory has been 
developed into an extremely powerful tool that has touched many areas of mathematics, 
science and engineering. Consider the following general problem of optimal control. 

2.1.1. The Optimal Control Problem 

Determine the state-control function pair, t h-)■ (x, u) G x that minimizes 
the cost functional 

J[x{-),u{-)]= j F{x{t),u(t))dt +E{x{tf)), (2.5) 

Jto 

subject to the dynamics. 


= f{x{t),u(t)), 


( 2 . 6 ) 
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initial conditions, 


x{to) = Xo, 


(2.7) 


at specified time, to, and endpoint conditions, 

e{x{tf)) = 0 , ( 2 . 8 ) 

where the running (or Lagrange) cost F : x —)■ M, the endpoint (or Mayer) cost, 

E : X R^^ —)■ M, / : R^^ x —)■ R^^ and e : R^^ x R^^ —)■ are Lipschitz 

continuous with respect to their arguments. A set of necessary conditions must be met in 
order to find candidate solutions to problem (2.5)-(2.8). Pontryagin’s Minimum Principle 
provides the necessary framework. 

2.1.2. Pontryagin’s Minimum Principle 

Pontryagin’s Minimum Principle was proved by Pontryagin in 1956 [22, 23]. It 
provides conditions that must be met in order for a solution to be considered optimal. As 
with the calculus of variations both necessary and sufficient conditions for optimality may 
be established. Although sufficient conditions are beyond the scope of this dissertation 
(see [24]), first order necessary conditions will be outlined with help from the calculus of 
variations. 

2.I.2.I. Calculus of Variations 

In the calculus of variations, problems of the form (2.1)-(2.2) are solved by consid¬ 
ering the variation of J, or A J, given by 

^J[y\y] = J[y]-J[y*] 

where y* is the minimizing curve, and y are all other admissible curves. For y* to be a 
minimizing curve it is necessary that /S.J[y*,y] > 0. Additionally, if all the first order 
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terms are collected in the expansion of A J, it is necessary that this collection (called the 
first variation or shown symbolically as 5J[y*]) must be equal to zero [21]. This same 
approach may be used to define the optimality conditions for problem (2.5)-(2.8). 


2.I.2.2. Necessary Conditions 

In order to apply a variational approach to problem (2.5)-(2.8), consider the aug¬ 
mented functional, 

A(-), z^(-)] = [ {F{x{t),u{t)) + X^{f{x{t),u{t)) - x(t))) dt 
Jto 

+E{x{tf)) - iy^e{x{tf)), 

where \{t) G and u{t) G are Lagrange multipliers, and \{t) is typically given the 
name costate or adjoint covector. As in the calculus of variations approach, considering the 
first variation, 6Ja[u*] = 0, a set of necessary conditions can be obtained [21, 25-27] 

m = ( 2 . 9 ) 

AW = -^. (2.10) 

where the Hamiltonian, H, is given by 

H{x{t),u{t),\{t)) = F{x{t),u{t)) + \^ f{x{t),u{t)). (2.11) 


Additionally, the Hamiltonian (2.11) reaches its minimum with respect iou aiu = u*. This 
is called the Hamiltonian Minimization Condition and can be expressed as 

u* = arg max i7(a;(t), M(t), A(t)), (2.12) 

uSU 
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where U defines a region of feasible control. Finally, the following conditions must be 
satisfied on the boundary 



(2.13) 


(2.14) 

x(tf)) =0, 

(2.15) 


where the endpoint Lagrangian, E, is given by 


E{x{tf), u) = E{x{tf)) + iy'^e{x{tf)). (2.16) 

Equations (2.9), (2.10), (2.12) and (2.13)-(2.15) provide the first-order necessary condi¬ 
tions for optimality and create the framework for Pontryagin’s Minimum Principle. 

2.I.2.3. Pontryagin’s Minimum Principle 

Lemma 2.1 (Pontryagin’s Minimum Principle). [26] Let, {x*{t),u*{t)), be a solution to 
problem (2.5)-(2.8). Then in order for x*(t) and u*(t) to be optimal, it is necessary that 
there exists a costate, \, and covector, o, that satisfies conditions (2.9), (2.10), (2.12) and 
(2.13)-(2.15). 

Remark 2.1. For problem (2.5)-(2.8), with added path condition, the following mixed 
state-control inequality path constraint is included, 

h(x(f),u(t)) < 0, (2.17) 

where h : x —)■ is Lipschitz continuous with respect to x and u. 
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With the addition of the path constraint, candidate solutions to problem (2.5)-(2.8) 


and (2.17) can be found by solving the nonlinear programing (NLP) problem 

u* = axg max H(x(f),u(f),\(f), p(f)), 

u&](x) 

where the constraint set, U C is given by 

I[J(a;) = {m G l]\h{x{t),u{t)) <0,x E G [to,tf]}. 

The augmented Hamiltonian (or Lagrangian of the Hamiltonian), H, is given by 

H(x(t),u(t),\(t),p(t)) = H(x(t),u(t),\(t)) + p^h(x(t),u(t)), 

where p(f) G are Lagrange multipliers. The modified set of necessary conditions 
are [27, 28] 


, dH 

w N dH 


along with the following conditions on the boundary 


K'^f) 

Hitf) 


dE 


dxf 


dE 

dtf 


e(x(tf)) = 0 . 


Additionally, the complementary (slackness) condition. 


(2.18) 

(2.19) 


( 2 . 20 ) 

( 2 . 21 ) 

( 2 . 22 ) 
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(2.23) 


<0, hi{x{t),u{t)) = 0, 

< 

1^=0, hi{x{t),u{t)) < 0, 

must be satisfied. 

Equations (2.18)-(2.23) provide the first-order necessary conditions for optimality 
for problem (2.5)-(2.8) and (2.17). 

Although Pontryagin provided a framework for finding candidate optimal solutions, 
many problems of optimal control are too difficult to solve analytically. It is easy to see the 
difficulty in solving the 2Nx Hamiltonian system of differential equations (2.9)-(2.10) or 
(2.18 )-(2.19). For this reason, numerically methods have become extremely important in 
solving optimal control problems. 

2.1.3. Numerical Methods for Optimal Control 

Many numerical techniques have been investigated for solving optimal control prob¬ 
lems since Pontryagin proved the Minimum Principle in 1956. These optimal control meth¬ 
ods take two main forms, indirect and direct. Recent surveys of these techniques are pro¬ 
vided by Betts [29, 30], Trelat [31] and Ross [32] and a historical perspective by Stryk et 
al. [33]. Indirect methods (such as the shooting and multiple shooting methods) solve Pon- 
tryagin’s necessary conditions for optimality. Although these methods have been shown 
to solve a wide range of problems with great accuracy, they prove to be difficult to im¬ 
plement, due to the knowledge of the calculus of variations required and the difficulty of 
providing good initial guesses. In contrast, the direct methods (such as Euler, Runge-Kutta 
and collocation methods) discretize the cost function, problem dynamics, etc, at specified 
time points. Due to the fact that direct methods require no knowledge of the necessary con¬ 
ditions for optimality, and the accuracies that may be obtained, they have recently gained 
much attention. Of the direct methods, specifically the global orthogonal collocation meth¬ 
ods (a.k.a. pseudospectral methods) have proven to solve difficult problems with great ac- 
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curacy [4, 8, 34, 35] after becoming an actively researched topic in the 1990s by Elnagar et 
al. [1] and Fahroo et al. [2]. 

Pseudospectral (PS) methods for optimal control discretize the problem at speci¬ 
fied nodes, called collocation points. Due to the properties of the orthogonal family of 
collocation points (such as those found via the Legendre or Chebyshev polynomial ba¬ 
sis) approximations converge at spectral rates [6]. The most widely used Legendre PS 
method [1, 5, 6, 9, 10] is based on the LGL points [36]. However, the Legendre PS method 
may be based upon LGR or LG nodes as well [36, 37]. PS methods for optimal control 
have been formally implemented in the MATLAB-based software package DIDO [38] and 
NASA’s Fortran-based package OTIS [39]. 

There are four parts to the numerical solution to an optimal control problem us¬ 
ing a PS method: discretization of the system dynamics, discretization of the state-control 
constraints, integration of the cost function and solving the nonlinear program (NLP). The 
mathematical background associated with the first three steps will be discussed in the fol¬ 
lowing sections. Spectral methods are attractive for discretizing the problem’s dynamics 
due to their superior accuracy. Two global spectral methods, collocation and Galerkin, 
will be outlined in Section 2.3.1. Additionally, Galerkin methods may be formulated as 
element-based methods. These local spectral element methods will be outlined in Sec¬ 
tion 2.3.2. A fundamental task in the formulation of these global and local methods is the 
selection of good discretization points and the use of interpolating functions. Both will be 
discussed in detail in Section 2.2. Finally, numerical integration, or quadrature, is typically 
used to integrate the cost function and will also be outlined in Section 2.2. 

The resulting NLP can be solved by using a commercial sequential quadratic pro¬ 
gramming (SQP) software packages such as dense NLP solver NPSOL [40] and sparse NLP 
solvers SNOPT [41, 42] and SPRNLP [43]. A feasible solution can be found that satisfies 
the tolerances specified in the optimization problem by adjusting the order of polynomial 
used in the approximation. 
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2.2. Interpolation and Numerical Integration 

Interpolation and numerical integration serve an important role in the methods out¬ 
lined in this dissertation and will be discussed in this section. The general structure of this 
section follows that provided by Giraldo in Chapters 4 and 5 of [44]. 

2.2.1. Interpolation 

Polynomial interpolation is the method used to construct an iV-th order polynomial, 
or interpolant, x^{t), that approximates a function, x{t). This is typically done by ensuring 
the interpolant passes through the iV +1 known points, {(fi, thata;(fi) = x^{ti), 

fori = 0,1,..., N. This may be accomplished by using a finite sum such as 

N 

!"(«) = ^4,(()*,, (2.24) 

j=0 

where {xj}^^Q are the expansion coefficients and are the basis functions. Defining 

the basis functions, modes (such as Legendre polynomials) leads to modal 

type of interpolation. However, defining the basis functions in a nodal fashion such that 
^jiU) = Sij, for i,j = 0,1,...,N, where 

[o, i^j, 

(such as Lagrange polynomials) leads to nodal interpolation. 

2.2.1.1. Modal Interpolation 

In modal interpolation, the basis functions, in Equation (2.24) are typi¬ 

cally orthogonal polynomials and the eigenfunctions of the singular Strurm-Liouville prob¬ 
lem. Commonly used polynomials are: Legendre, Chebyshev, Fourier and Jacobi. For this 
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discussion the Legendre polynomial, L{t), will be the foeus, defined by [45] 


L,(t) 


(-ly di 

2ij ! dP 




(2.25) 


and therefore will be the ehosen basis. The Legendre polynomials result from the speeial 
ease of the singular Strurm-Liouville problem [45], 


Figure 1 shows the first seven Legendre polynomials. The speetral eoeffieients, 



Figure 1: Legendre polynomials, 


for the eontinuous Legendre expansion, are defined as [46] 


ttj = — / x{t)LAt)dt, 

Aj J-i 


(2.26) 
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where the normalizing constants, {7j}^0’ the Legendre polynomials are given by 

2 

~ 2j + l 

The truncated Legendre modal expansion is then given by 

N 

j=o 

However, due to the interpolatory nature of x^, it is natural to seek spectral coefficients, 
defined by 

N N 

(2.29) 

j=0 j=0 


for the known points where V is the generalized Vandermonde matrix given by [47] 



^ Lo(to) 

Li{to) 

LAr(to) ^ 



V = 

Lo{ti) 

Li{ti) 

Lfq(ti) 


(2.30) 


y Loitn) 

Liitn) ■ 

■ ■ L]s[(tN) y 



From Equation (2.29), the modes, {aj 

}^ 0 ’ ^i^d the nodes, 

are related by the 

generalized Vandermonde matrix (2.30) by the relationships [47] 




N 

= Yy^aj 

i=o 

and 

N 

( 2 . 31 ) 

j=0 


(2.27) 


(2.28) 
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where = x^{tj) for j = 0,1,..., iV. 

Remark 2.2. Note that to form Equation (2.31) the Vandermonde matrix (2.30) must be 
invertible and therefore nonsingular (and more practically speaking, well-conditioned). We 
will see that the invertibility of the Vandermonde matrix is dependent upon the interpolation 
quality of grid, Section 2.2.1.4 for grid selection) [47], Additionally, we take 

comfort in the fact that the set of Legendre polynomials, orthogonal system 

that has shown to produce well-conditioned Vandermonde matrices for carefully selected 
nodes (as compared with the ill-conditioned Vandermonde matrices of non-orthogonal sys¬ 
tems such as the the power basis, [48]. 

Note that Equation (2.28) is a sum of frequencies, and amplitudes, 

that together compose the (N-yl) modes of x^ (t). It is thus fitting to describe this approach 
as modal interpolation. 


2.2.I.2. Nodal Interpolation 

In nodal interpolation, the basis functions, iii Equation (2.24) are the Ea- 

grange polynomials, of order N, defined on grid obtained from the gen¬ 

eral definition [45] 

< 2 - 32 ) 

i=o Vj V 
i4j 


Additionally, the Eagrange polynomials may be defined in terms of the Eegendre polyno¬ 
mial by [46] 



1 (t^-l)Ljv(t) 
iV(iV + l) (t-tfL^(tf 


(2.33) 


Eigure 2 shows the order N = 6 Eagrange polynomials, defined on an equi- 

spaced grid, f e [—1,1]- 
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Figure 2: Lagrange polynomials of order N = Q defined on an equi-spaeed grid. 

The nodal interpolation of the function x{t) can be accomplished by the iV-th order 
expansion 

N 

x^it) = (2-34) 

j=0 

where = x^{tj), for j = 0,1,..., iV, since the Lagrange polynomial has the property, 

(L) = 5ij. 

Additionally, the Legendre polynomial, Li{t), of order i can be written as linear 
combinations of Lagrange polynomials, of order N defined on grid, {UjfLo, 

by the relationship [47] 

N N 

U(t) = (2.35) 

J=0 j=0 
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Likewise, the Lagrange polynomial may be written as linear combinations of Legendre 
polynomials by the relationship 


N N 

= E hit) = (2.36) 

j=0 j=0 

From Equations (2.31) and (2.36) we can relate the modal and nodal forms of the interplant, 

x^{t), by 


N 

i=0 

N / N 


N / N 


i=0 \j=0 


N 


X 


Ni 


=Y,h(t)[J2hh:t:f‘]=Y,h(t)a,. 

j=0 V *=0 / j=0 


Therefore, Equation (2.34) is truly a nodal representation of Equation (2.28). 


2.2.I.3. Transformations between grids 

Consider the problem of transforming between two different grids, and 

where M < N. Eet the set of Eagrange polynomials of order M 

defined on grid {rjjjiQ. Also, let the function x^{t) be the Eagrange interpolating poly¬ 
nomial 

M 

= E'^f 

, 7-0 

where x^^ = x^{Tj), for j = 0,1,..., A^, since 0^(ri) = 6ij. Then the approximation of 
x^ at the dense gridpoints, {tk}k=o, he calculated with the linear transformation 

M M 

! = 0, 1, . . . , AT, 

j=0 j=0 
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where the (iV + 1) x (M + 1) linear mapping, is given by 


/ 


rjnNM _ 




0f(to) 


mo) 

mi) 


\ 


(2.37) 


\W{tN) 01 (^Ar) ••• 4>M{tN) J 


In a similar fashion, the approximation of x may be transformed between two grids. 
Note that the approximation of x on grid {rj }^o be given by 

M 

x{t )«i"(«) = 

where the derivative of the Lagrange polynomial is defined as 


M 


?'«) = E 


3 


k=0 


M 

n 


t — tj 


k¥=j 


i^k 


(2.38) 


Then the approximation of x^ at the dense gridpoints, can be calculated with the 

with the linear transformation 


M M 

j=0 j=0 


where the (iV + 1) x (M + 1) linear mapping, , is given by 


^ 0^(^o) 0f(fo) ••• 0M(fo) ^ 

05^(ti) 0f(ti) mi) 
V0O^M 0rM ■■■ 0MMy 


(2.39) 
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Using the transformation matrices (2.37) and (2.39) to relate different grids will 
serve as an important tool for the multi-scale approximation methods outlined in Chapters 3 
and 5. However, the accuracy of interpolation is extremely important and will be discussed 
next. 

2.2.I.4. Interpolation Quality 

Approximation quality is a great concern when using interpolation. The goodness 
of the approximation x^{t) is directly related to the grid points, from which the 

Lagrange polynomials, are defined. A measure of interpolation goodness is the 

Lebesque constant, Aat, given by [45] 

N 

Ajv = maxy^Uf(t)|. (2.40) 

The best interpolating polynomial x^{t) is one that minimizes the Lebesgue constant (2.40), 
due to the following result [45], 

\\x{t)-x^{t)\\^^ < (1 + A 7 v)||a;(t) -p(f)||^oo, 

where p{t) is the best approximating polynomial of x{t) in the L°°-norm (see Appendix A). 
From [45, 49], for any set of {N + 1) distinct points, U G [—1,1], for f = 0,1,..., A^, the 
Lesbegue constant (2.40) has the lower bound [45], 

2 

-log(A^ + 1) + a < Aat, 

TT 

where a = f (7 + log^) ~ 0.521 and 7 = 0.57721566... is the Euler-Mascheroni con¬ 
stant. So, at best the selected grid is associated with a Lebesgue constant that grows 
logarithmically [45]. Common Legendre family of points used for interpolation are the 
Legendre-Gauss (LG), Legendre-Gauss-Lobatto (LGL) and Legendre-Gauss-Radau (LGR) 
points. 
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Legendre-Gauss Points. The LG points, are defined by —1 < to < ■ ■ ■ < < 1, 

and are the roots of 

e(t) = L^+i(t), (2.41) 

where LAr+i(t) is the {N + l)-th order Legendre polynomial. Note that the LG points do 
not inelude the endpoints, t = ±1. Figure 3 shows the LG points for various orders of N. 


\ .^a|/\|/ s|/ \|/ 

X 1 


- * - 

(- ^ 

^ 

^ 1 


- 

^ ^ 

-^ 

— 

$- 1 

* - 

*- 

— 


Figure 3: LG points for N = 10, 20 and 30. 


Legendre-Gauss-Lobatto Points. The LGL points, {tj^o’ defined by to = —1 < 
ti < • • ■ < tAT = 1, and are the roots of 

e(t) = (l-t2)LAr(t), (2.42) 

where L^it) is the derivative of the N-th order Legendre polynomial. Note that the LGL 
points inelude the endpoints, t = ±1. Figure 4 shows the LG points for various orders of 
N. 
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Figure 4: LGL points for N = 10, 20 and 30. 


Additionally, Figure 5 shows the order N = 6 Lagrange polynomials, 0^, defined on a 
LGL grid. 



Figure 5: Lagrange polynomials of order N = Q defined on a LGL grid. 

Legendre-Gauss-Radau Points. The LGR points, are defined by Iq = —1 < ti < 

■ ■ ■ < tpf < 1, and are the roots of 

+ L^it). (2.43) 
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Note that the LGR points only inelude the endpoint, t = —1. Figure 6 shows the LGR 
points for various orders of N. 




N <!✓ xi/ \l^ nL' xI/ nU nI^ \1x xI/ \l' \l' xl/ xl^ xly xLxL • 

Tfr” "Tfr - "Tfr "Tfr ■ 7]^ 7^ - Tfr 7^ - Tfr Tfr -Tfr Tfr yf^ Tfr TfrTfTf. 




-0.5 0 Xs 1 

Figure 6: LGR points for N = 10, 20 and 30. 


^ 


Additionally, flipped-LGR (F-LGR) points are the negative of the LGR points and are there¬ 
fore defined by—1 < to < ■■■ < = 1- Note that the F-LGR points only inelude the 

endpoint, f = 1. 

Although all three sets of Legendre points (LG, LGL and LGR) have Lebesgue 
constants (2.40) that grow logarithmically or sublinearly with N, the LGL grid is asymp¬ 
totically associated with the near optimal Lebesgue constant [45], 


< -\og{N + 1) + 0.685.. 

TT 


As alluded to in Remark 2.2, the quality of interpolation can also be observed by 
analyzing the conditioning of the Vandermonde matrix (2.30). Due to the relationship 
between the Legendre polynomials, {Lj}^Q, and Lagrange polynomials, {(j)^}f^o, shown 
in (2.35), Cramer’s rule [50] provides the following relationship 


^ _ Det [L(fo),..., L(tj_i), L(f), L(fj+i),..., L(f7v)] 
^ Det [V^] 


(2.44) 


where L(f) = [Lo(f), Li(f), ..., L 7 v(f)]^. As pointed out by Hesthaven et al. [47], if the 
goal is to minimize the Lebesque constant (2.40), we should strive to maximize the de¬ 
nominator of Equation (2.44), Det[L^]. This leads to the LGL grid set [51]. Additionally, 
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the Chebyshev-Gauss family of points proves to have excellent interpolation quality, par¬ 
ticularly the Chebyshev-Gauss-Lobatto points when measured by the Lebesgue constant 
growth [45, 52]. However, the focus in this dissertation will be on the Legendre basis, and 
thus the LG, LGL and LGR points. 

Unfortunately, equi-spaced points prove to be a very poor grid selection for inter¬ 
polation. The Lebesgue constant for the equi-spaced points grows asymptotically like [52, 
53], 

oTV-l-l 

aES ^ ^ _ 

^ eN(\og iV 4 - 7 ) ’ 

very far from optimal. 

As an example of interpolation quality consider the function 

f(t) = cos(/i 7 rf), t e [-1,1], (2.45) 

with /i = 3, shown in Figure 7. 



Figure 7: Plot of f(t) = cos( 37 rt). 
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When low order approximations are used to interpolate f(t), inaeeuraeies are appar¬ 
ent. Figure 8 shows the inaeeuraeies in the 10-th order Lagrange interpolating polynomial 
approximation of f{t) with LGL points. 



Figure 8: Interpolation of f{t) = cos^Snt) with 10-th order LGL points. 
However, as the interpolation order, N, is inereased, the maximum error, 

||error||^ = ||/(^i) “ IL’ * = 0,1,..., iV, 

deereases exponentially with N, where || ( ||oo represents the maximum element of veetor, 
G M". Figure 9 shows the visual aeeuraey of the 30-th order Lagrange interpolating 
polynomial of f{t) with LGL points. 
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1.5 


— Exact 

— Numerical 



-1.5'-^^-'- 

-1 -0.5 0 0.5 1 

t 


Figure 9: Interpolation of f(t) = cos(37rt) with 30-th order LGL points. 


Additionally, Figure 10 eompares interpolation of f(t) with equi-spaeed, LG, LGL 
and LGR points, for various orders of N. Notiee that for LG, LGL and LGR points, the 
maximum interpolation error drops to 0(10“^®) by iV = 35. However, in general, the 
equi-spaeed points prove to have very poor interpolation quality. 
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Figure 10: Comparison of interpolation errors for various orders of N and equi-spaced, 
LG, LGL and LGR points. 


Due to the accuracy of LGL interpolation, and inclusion of the endpoints, t = ±1, 
LGL points are readily used for numerical computation. However, also related to point 
selection is the accuracy of numerical integration. This is an important factor for the direct 
methods for optimal control (due to the cost function that is normally approximated by 
numerical integration) and will be discussed next. 


2.2.2. Numerical Integration 

Numerical integration, or quadrature, is a way of approximating an integral with a 

sum 



N 

x{t)dt ^ y^x{tk)wk, 

k=0 
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where {wk}^^Q are the quadrature weights and {tk}k=o the associated points. Ideally, 
the numerical integration is exact, but it is reasonable to expect an error, e^, such that 

„1 N 

= / x{t)dt - y^x{tk)wk. 

k=0 

For the special case of the numerical integration of general function x G that is 

approximated by the Lagrange interpolating polynomial 

N 

x{t) ^ X^{t) = '^(j)f{t)x^\ 
j=0 

the error may be given by [54] 


e^ = 


iN + l)\J_ 


„1 N 

/ Ylit-tj 

"'-1 j =0 


' di(^) 


dt, 


for arbitrary function ,^(t) e [—1,1]. 

However, for certain classes of polynomial functions, the quadrature error is zero. 
Consider the function x{t) G Pn, represented as the finite sum 


N 

x{t) = 

j=0 

Performing numerical integration on x results in 




{t)dt 


N N N 

= {tk)wk = y^x^^wj. 

j=0 fc=0 j=0 
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Clearly, since 



(2.46) 


quadrature is exact \/x{t) G Pn- Additionally, the quadrature weights, can be 

found with the relationship 


Wj = 



(2.47) 


2.2.2.1. Gaussian Quadrature 

Consider now the case that x{t) G P 2 N +1 written in the form 


x{t) = LN+i{t)f{t)+g{t), 


where f,g E Pn, Ln+i is the Legendre polynomial of order {N + 1) and Pn denotes 
the space of all polynomials of degree < N. Also consider the {N + 1) points, 
that are the roots of LN+i{t) (known as LG points, discussed in Section 2.2.1), and the 
associated quadrature weights, (known as LG quadrature weights, found via the 

general definition (2.47) or the more specific definition (2.48)). Then x{tk) = g{tk), for all 
k = 0,1,..., N, and 

j x{t)dt = j {LN+i{t)f{t) + g{t))dt 
N N 

= J2LN+i{tk)f{tk)wk + J2g{t k)Wk 
k=0 k=0 

N N 

= y^,g{tk)wk = 'y^x{tk)wk. 

k=0 k=0 

This is known as LG quadrature (or simply Gauss quadrature), which is exact Wx{t) G 
P 2 N+ 1 - In the case that the function x{t) is a polynomial, such that a;(t) G P 2 N+ 5 , numerical 
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integration is exact for LGL and LGR quadrature, where <5 = — 1 and 0, respectively. The 
proof of LGL and LGR quadrature exactness for polynomials is similar to that given above 
for LG quadrature. The list of LG, LGL and LGR weights are presented below (provided 
by [47]) and point locations are given by Equations (2.41), (2.42) and (2.43), respectively. 
Legendre-Gauss Quadrature. The LG quadrature weights, are given by 


[l-(4)2][L^+i(4)p- 

Legendre-Gauss-Lobatto Quadrature. The LGL quadrature weights, {wkj^^Q, are given 
by 


Wk 


2 1 

iV(iV+l)|LA,(4)l"' 


(2.49) 


Legendre-Gauss-Radau Quadrature. The LGR quadrature weights, {wk}k=o^ given 
by 


{N + 1)2 [L7v(tfc)]^ 


(2.50) 


Additionally, the F-LGR quadrature weights, {wk}k=o, ^be reordered LGR quadrature 
weights, {wk}k=o = {wN-k}k=o- 


2.2.2.2. Gaussian Quadrature Accuracy 

The Legendre-Gaussian family of quadrature is widely used due to its integrating 
accuracy and has the following error estimates for LG, LGL and LGR quadrature, respec¬ 
tively [55] 
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^LG 


^(2iV + 3)[(2iV + 2)!]3 cit(2^+2) 


, ee [-1,1], 


^LGL 


-(TV + l)iV322^+i[(iV - 1)!]^ 


(2iV + l)[(2iV)!]3 




, ee [-1,1], 


cN 

-LGR 


22iv+i(jY ^ i)(iV!)4 
' ;(2iV +1)!]3 cit(2^+i) 


, ee [-1,1]. 


Consider again the function (2.45), with /i = 1/2. As an example of quadrature 
accuracy, consider the numerical approximation. 


cos 


'-1 


Tit 


N 


dt ^ y^j^{tk)wk, 


(2.51) 


A:=0 


where 


N 

/(i) =» /«(«) = 


j=0 


and = cos(^), for j = 0,1,..., A^. Due to the property of the Lagrange polynomial, 
= 6ij, the quadrature expression in (2.51) takes the form (2.46), 


N 


N 


5^/^(4)wfc = 


Wk. 


k=0 


k=0 


As with interpolation error, quadrature error decreases exponentially with N. Fig¬ 
ure 11 shows numerical integration error. 


error = 


-1 


cos(y )dt - y^cos{^)wk 

1 k=0 
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for LG, LGL and LGR quadrature, all of which demonstrate similar exponential conver¬ 
gence as N increases. 



Figure 11: Comparison of quadrature errors for various orders of N and LG, LGL and LGR 
points. 


Due to the accuracy of quadrature, quality of interpolation and inclusion of the 
endpoints, t = ±1, LGL points are an important part of many numerical computation ap¬ 
plications, to include direct methods for optimal control. For instance, in the Legendre PS 
method the problem state-control constraints are discretized using LGL points. Addition¬ 
ally, the states are approximated with globally interpolating Lagrange polynomials defined 
on an LGL grid and the cost function is approximated using LGL quadrature rule. LGL 
interpolation and quadrature also serve an important role in the construction of the Galerkin 
optimal control formulations discussed in this dissertation. 

A crucial step in solving optimal control problems with direct methods—and yet to 
be discussed—is the discretization of the system dynamics. Although a number of tech¬ 
niques have been investigated to do this, spectral methods have proven to be effective and 
efficient. This is a highlight in the PS direct methods for optimal control, where collocation 
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methods are used to discretize system dynamics. Spectral methods also serve as the heart 
of Galerkin optimal control, where Galerkin methods are employed. Both collocation and 
Galerkin methods will be discussed in the next section. 

2.3. Numerical Solutions to Differential Equations 

Spectral methods, which have gained much popularity due to their spectral accu¬ 
racy and versatility [56, 57], can be formulated for both local (element-based) and global 
approximations. Consider the task of discretizing the dynamics of problem (2.5)-(2.8), 

x{t) = f{x{t),u{t)), t e [to,tf], (2.52) 

where t h-)■ {x, u) G x with initial conditions. 


x{to) = Xo, 


(2.53) 


and endpoint conditions. 


e(a^(f/)) = 0, (2.54) 

where / : x —)■ and e ; x —)■ are Lipschitz continuous with 

respect to their argument. This can be accomplished using a number of spectral method for¬ 
mulations. However, two methods will be discussed here, collocation and Galerkin. Addi¬ 
tionally, of the spectral element methods, continuous and discontinuous Galerkin element- 
based formulations will be outlined. 
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2.3.1. Spectral Methods 

The starting point for the spectral approximation of Equation (2.52) is to approxi¬ 
mate the solutions x and u by the finite sums 

N 

j=0 

N 

j=0 

where Xj and Uj are expansion coefficients and and are basis functions. In terms of 
the approximation x^ and u^, Equation (2.52) becomes 

*"«) - «"(?)) = e"K). (2.55) 

where At = tf — to and is the error (or residual) in the approximation which, gener¬ 
ally, is not zero. The relationship between the physical time domain, t G [to,tf], and the 
computational space, ^ G [—1,1], is given by 

2 2 
f = Y— (t — to) — I and df = ^—dt, 

^ At ^ ^ At 

and conversely. 

At At 

t = — {^ + 1) + to and dt = -^d^- 

2.3.1.1. Collocation 

In the collocation method the basis functions, {'I’jjjlo the Eagrange polyno¬ 
mials (2.32), {0^}^Q, of order N, defined on the grid of collocation points, ^ 

[—1,1]; while where is any continuous function (not nec¬ 
essarily a polynomial) with the property = dij. The expansion coefficients are, 

Xj = x^^ and Uj = , therefore x^{Q) = x^^ and u^{Q) = for j = 0,1,..., iV. 
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The approximations of x and u in the computational space, are given by 

N 

j=0 

N 

j=0 

The approximation of x in the computational space, is then 

N 

j=0 


where the derivative of the Lagrange polynomial, given by Equation (2.38). 

Additionally, in the collocation method, the error term, e^, is ideally forced to zero at each 
collocation point, therefore. Equation (2.55) becomes 


/(s'". fi“)=o. 

j=0 


i = 0,1,..., N. 


Collocation methods assume the title pseudospectral methods, due to the nodal 
nature of the formulation—in lieu of spectral referring to a transformation from physical 
to spectral space. Common Eegendre family of collocation nodes used are the EGE, EG 
and EGR points. When EGE nodes, defined by —1 = = 1, 

are used for the discretization, the Eegendre PS differentiation matrix. A, is given by Aij = 
(j)^ (Ci) for hj = 0,1,..., N, resulting in the system of equations. 


^ .At . . 

'^AijX^^ - —f{x^\ =0, i = 0,1,..., A^. (2.56) 

i=o 
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In addition to Equation (2.38), the Legendre PS differentiation matrix, A, may be defined 
in terms of the Legendre polynomials by [46] 


Aij 


= < 


LN((i) 1 
LnH]) ii-ij 

N{N+1) 

4 

Af(7V+l) 

4 

0 


i ^ J, 
i = j = 0, 
i=j = N, 


(2.57) 


Remark 2.3. The derivative of at each LGL point exactly equal to 


N 


'(f.) = 

j=0 


for any polynomial with degree less than or equal to N [46]. However, a feasible solution 
to the equality dynamical constraint may not exist. In order to guarantee feasibility of the 
discretized problem, Gong et al. [5], suggest a relaxation of the equality constraint. 


Therefore, Equation (2.52) may be diseretized with the following inequality eonstraint, 

where 5^ is the feasibility toleranee that is dependent on N and the smoothness of x and 
u (see Seetion 3.2); and || C ||oo represents the maximum element of veetor, ( G MA. The 
initial eonditions and endpoint eonditions may be approximated similarly by 




\x 


7V0 


— xn <5^ and e(a; 

^ rv~) — ^ 


7X.NN 


< 5 


N 


Colloeation methods have beeome popular for the diseretization of system dynam- 
ies in direet methods for optimal eontrol, speeifieally psuedospeetral methods. Lor the 
Legendre PS method, the LGL points beeome the diseretization of ehoiee. 
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2.3.I.2. Galerkin Numerical Methods 


Galerkin methods can be subdivided into two main categories, Bubnov-Galerkin 
and Petrov-Galerkin. The weighted residual method forms the basis for the development 
of these approximation techniques and will help to distinguish them [57]. 

For the weighted residual method, the weak integral form of the Equation (2.55) is 
solved by multiplying by a test function, integrating over the domain, and ideally we 
force the residual term to zero, 

«"(?))) d( = = 0, (2.58) 

for i = 0,1,..., iV. 

Remark 2.4. Setting the residual terms to zero in Equation (2.58), 

for each i = 0,1,..., N), is akin to forcing the orthogonality of the space spanned by '1/j 
and in L^[—1,1]. 

The approximation and can be found satisfying Equation (2.58). Common 
test functions include orthogonal polynomials (such as the Eegendre polynomials, L) and 
the trigonometric functions. In the Bubnov-Galerkin method, the test functions are the 
same as the basis functions, unlike the Petrov-Galerkin method, where the test and basis 
functions are different. The general structure of these global Galerkin methods—as well 
as the mathematical notation used in this dissertation—is provided by Giraldo [44] and 
discussed in the following sections. 

Bubnov-Galerkin In the Bubnov-Galerkin method (or often called simply the Galerkin 
method) the test and basis functions are the same. These functions can be modal or nodal 
in nature, however, in this section the focus will be on nodal Galerkin methods. Eor a nodal 
Galerkin approach, it is common to use a Eegendre based grid such as the EGE, EG or EGR 
nodes, and define the test and basis functions as Eagrange polynomials (2.32), 
of order N, on the selected grid. A popular selection for interpolation points are the EGE 
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nodes due to the aeeuraey of LGL quadrature and the inelusion of endpoints, t = ±1. For 
this diseussion, the LGL nodes will be the foeus. 

The approximations of x and u in the eomputational spaee, ^ G [—1,1], are given 
by 

N 

3=0 

N 

3=0 

where is any eontinuous funetion (not neeessarily a polynomial) with the property 

= Sij. The expansion eoeffieients are, xj = x^^ and Uj = , therefore x^{^j) = 

x^^ and for j = 0,1,..., iV. Equation (2.58) beeomes [44] 

J 

for i = 0,1,..., iV, or using matrix-veetor notation, 

N 

'^DijX^^ -Ci = 0, i = 0,1,... ,N. 

j=0 

The Galerkin differentiation matrix, D, and RHS veetor, c are defined as 

A, = r A" (?)<(?)<*?. 

fori,j = 

Using LGL quadrature, D can be ealeulated with the relationship 

Q 

dd>i3 = = AjWi, i,j = 0,1,... ,N, 

k=0 
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where are the LGL weights given by Equation (2.49) and A is the Legendre PS 

differentiation matrix (2.57). Since LGL quadrature rule is exact for polynomial integrands 
of degree less than or equal to 2N — 1, the numerical integration is done exactly when 
Q = N LGL integration points are used. 

Using LGL quadrature rule, c can be approximated by the relationship 

Af ^ 

Ci ~ i = 0,1,..., iV. 

k=0 

When Q = N LGL quadrature points are used, the RHS vector approximation, c^, can be 
expressed in the simplified form 

= i = 0, 

Remark 2.5. Recall that for LGL quadrature rule, integration is exact for polynomial 
integrands of degree less than or equal to 2N — 1. IfQ = (iV + 1) integration points are 
used, the RHS vector will integrate exactly when f{x(t),u(t)) is linear in x and u. In the 
case of a nonlinear function f, the accuracy of integration (and therefore the accuracy of 
the overall approximation) can be improved by increasing the number of quadrature points 

Q. 

When Q = N LGL quadrature points are used to calculate the Galerkin differentiation 
matrix and approximate the RHS vector, the system may be simplified as 

N 

= 0, i = 0,1,..., iV. (2.60) 

j=0 

Remark 2.6. Inform (2.60), the resulting Galerkin equations that must be satisfied are 
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for i = 0,1,..., iV. Note that the relationship in parentheses, 


j=0 

for i = 0,1,..., N, are the same equations that would be satisfied when using the col¬ 
location method (see Equation [2.56]). For this reason, Bubnov Galerkin with numerical 
integration is sometimes called the “collocation method in the weak form.” [57] 

In the words of John Boyd, “collocation—with the right set of points—must inherit 
the aura of invincibility of the Galerkin method.” [58] 

Remark 2.7. An inequality version of (2.61) has been known and used in pseudospectral 
optimal control methods. Details on its relationship with Galerkin optimal control are 
addressed in Chapter 4 in Remark 4.2. 

Due to the results of Gong et al. [5], we know a feasible solution to the equality 
dynamieal eonstraint may not exist. In order to guarantee feasibility of the diseretized 
problem, the following inequality eonstraint is suggested, 

<5^, i = D,l,...,N, (2.62) 

where is the feasibility toleranee that is dependent on N and the smoothness of x and 
u (see Chapter 4). The initial eonditions and endpoint eonditions may be approximated 
similarly by 


N 

.7=0 


— Xoll^ < and ||e(a;'^^)||^ < 5^. 

Remark 2.8. The inequality formulation (2.62) introduces some fundamental differences 
in numerical analysis. In the Galerkin approach, the error is measured by the Lf-norm. 
As a result, 5^ has a feasibility with a slightly relaxed bound, by a factor of y/wi (see 


38 



Section 4.2 for the general Galerkin optimal control computational strategies, particularly 
Equations (4.10) and (4.13)). 


Remark 2.9. With the Galerkin formulation outlined here, the initial conditions may be 
enforced in a weak sense. In other words, ICs may be imposed only up to the order of ac¬ 
curacy of the numerical approximation itself. Consider again Equation (2.59). Integration 
by parts on the first term results in the Galerkin weak form. 



In terms of the approximating polynomials (and introducing the true initial condition, 
a;^(—1) —)■ x(—l)) we have 

_ /■! Ay- yl 


for i = 0,1,..., N. Integration by parts, yet again, results in the Galerkin strong form. 


N / N 

j=o \j=o 

( ^ 

- ] -Ci = 0. 

\j=0 


(2.63) 


Equation (2.63) may be formulated for weak enforcement of ICs by letting x(—l) = xq and 
a;^(l) = x^^ . Additionally, when Q = N LGL quadrature points are used to calculate 
the Galerkin differentiation matrix and approximate the RHS vector, the system may be 
simplified as 


N 


'^DijX^^ + Ki- 
j=0 


0 , 


(2.64) 
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for each i = 0,1,..., iV, where 


I — Xq, i = 0, 

[O, 

The IC term k now provides a natural way to introduce initial conditions into the discretiza¬ 
tion of the dynamics. Again, in order to guarantee feasibility of the discretized problem, 
the following inequality constraint is suggested, 

<6^, i (2.65) 

OO 

where 5^ is the feasibility tolerance that is dependent on N and the smoothness of x and u 
(see Section 5.1). Finally, the endpoint conditions may be approximated similarly by 

Remark 2.10. In [11], the equation resulting from dividing Equation (2.64) by Wi is in¬ 
troduced for primal-only closure conditions. However, for feasibility the inequality version 
of this expression. Equation (2.65), must be used for computational purposes. It should 
be noted that if the equation in [11] is multiplied by Wi first, then relaxed as an inequality 
bounded by 5^, the resulting inequality would be in agreement with the feasibility of the 
Galerkin weak boundary formulation discussed in Section 5.1.2 (see Equation (5.6)). 

Petrov-Galerkin In the Petrov-Galerkin method the test and basis functions are different. 
As with the Bubnov-Galerkin method, these functions can be modal or nodal in nature. In 
this section the focus will be on selecting a modal test function and a nodal basis. This 
will create the framework that will be used in Chapter 4. For this formulation, the selected 
test functions will be the Legendre polynomials, {Lj}^^Q, and the Lagrange polynomials 


N 

-4 Ki- 

7=0 
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(2.32), of order N, will be the basis. Again, for this discussion, the LGL node 

structure will be used for the problem discretization. 

The approximation of x and u in the computational space, G [—1,1], are given by 

N 

j=0 

N 

j=0 

where is any continuous function (not necessarily a polynomial) with the property 

= Sij. The expansion coefficients are, Xj = x^^ and uj = , therefore x^{^j) = 

x^^ and u^i^j) = for j = 0,1,..., A^. Equation (2.58) becomes [44] 

£ (0<ie S"" - Y f = 0. i = 0,1,..., Af, 

or using matrix-vector notation, 

N 

= z = 0,l,...,N. 

j=0 

The Galerkin differentiation matrix, D^, and RHS vector, are defined as 

At 

cf=^y 

forf,j = 0,1,...,A^. 

Using LGL quadrature, can be calculated with the relationship 

Q 

Ay = z,j = 0,l,...,N, 

k=0 
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where are the LGL weights given by Equation (2.49). Again, since LGL quadra¬ 

ture rule is exact for polynomial integrands of degree less than or equal to 2N — 1, the 
numerical integration is done exactly when Q = N LGL integration points are used. 

Using LGL quadrature rule, can be approximated by the relationship 

A/ Q 
fc =0 

When Q = N LGL quadrature points are used, the RHS vector approximation, c^, can be 
expressed in the simplified form 


^Ni 


At ^ 

A:=0 


Nk -Nk 


)wk, i = 0,l,...,N. 


Remark 2.11. Again, if Q = {N + 1) integration points are used, the RHS vector will 
integrate exactly when f{x(t),u(t)) is linear in x and u. If f is a nonlinear function, 
accuracy of integration may be improved by increasing the number of quadrature points Q. 

When Q = N LGL quadrature points are used to calculate and c^, the system 
may be simplified as 


N 


j=0 


0, i = 0,1,..., A^. 


( 2 . 66 ) 


A feasible solution to the equality dynamical constraint may not exist. In order to guarantee 
feasibility of the discretized problem, the following inequality constraint is suggested. 


N 

j=0 


-Ni 


<5^, i = 0,l,...,N, 


where 5^ is the feasibility tolerance that is dependent on N and the smoothness of x and 
u (see Chapter 6). The initial conditions and endpoint conditions may be approximated 


42 



similarly by 


||x^°-a;o|L <5^ and ||e(x^^) L < 5^. 

It is clear that in the Petrov-Galerkin formulation (2.66), the differentiation matrix, 
D^, and RHS vector, c^, do not simplify as cleanly as given for the Bubnov-Galerkin 
formulation (2.60). This inevitably will have negative effects on computational efficien¬ 
cies. However, casting the problem in the Petrov-Galerkin numerical form will have nice 
consequences when applied to Galerkin optimal control, as will be shown in Chapter 6. 

2.3.2. Spectral Element Methods 

Spectral element methods are local (elemental) applications of spectral methods. 
They combine the flexibility of finite elements with the accuracies associated with spectral 
methods. This element-based numerical approach is advantageous due to its ability to 
handle complicated geometries and can be easily formulated for adaptive strategies [16, 
59]. In this section, the focus will be on two Galerkin formulations, continuous Galerkin 
and discontinuous Galerkin element-based methods. Continuous Galerkin techniques were 
first applied to ordinary differential equations (ODEs) in 1972 by Hulme [14, 15] and a 
study of global error control was done by Estep et al. [60] in 1994. The first analysis of 
discontinuous Galerkin methods applied to ODEs was done in 1974 by Reed et al. [61] 
and an adaptive error control technique was used by Bottcher et al. [62] in 1997. More 
recently, multi-adaptive continuous Galerkin and discontinuous Galerkin techniques have 
been studied by Eogg and presented in a series of papers [63-65]. The general structure of 
these element-based Galerkin methods—as well as the mathematical notation used in this 
dissertation—is provided by Giraldo [44] and discussed in the following sections. 

2.3.2.1. Continuous Galerkin 

Consider a continuous element-based Galerkin approach to discretizing (2.52). Again, 
for this discussion, the EGE node structure will be used for the problem discretization. In 
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this approximation, the weak integral form of (2.52) inside each element, takes the 
form [44] 


[ dt = 0, (2.67) 

J S7e 

for e = 1, 2,..., iVe and i = 0,1, ..., N, where defines the total domain. 

The state trajectory, x{t), is approximated inside each element, fie, by interpolating iV-th 
order Lagrange polynomials, at the nodes by the relationship 

N 

j=0 

for e = 1, 2,..., iVe, where are the LGL nodes, mapped back to the 

physical space inside each element, fie- Also, let u^{t) be an interpolating function of 


N 


U 


(e)N 


(«) = 


i=o 


where (t) }^o are any set of continuous functions (not necessarily polynomials) with 


the property = Sij. Therefore for e = l,2,...,A'e and 


(e)\ 


j = 0,1,..., A^, and similarly, The relationship between the physical 


time domain, t G [to,tf] = 
by [44] 


AO ANe) 
I'O , I'M 


, and the computational space, ^ e [—1,1], is given 




At(e) ^0 


(e) 


1 and dp, = 




dt. 


and conversely. 




A) 

0 


and dt = —-— dp. 
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where — is the size of eaeh element, f2e, whieh ean be nonuniform in length. 

The Lagrange polynomial defined on the LGL eomputational domain is given by 


j¥=i 

The state trajeetory, x, can now be approximated inside eaeh element, fie, by 

N 

j=0 

where the Lagrange polynomials defined on the LGL grid. Likewise, u^{^) 

is given by 

N 

i=o 

where (6) = 

Remark 2.12. In this formulation = ^(e+i)Afo ^(e)NN _ ^{e+i)m^ ^ _ 

1, 2,..., Ae — 1. This continuity condition is a consequence of the global formulation of 
the problem discussed in Remark 2.13. 


In the eomputational domain, the system beeomes 


-1 Ay-(®) 


'-1 


'-1 


for e = 1, 2,..., Ae and i = 0,1,...,A, and in terms of the approximating polynomials 
beeomes 




j=0 


'-1 


'-1 
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In matrix-vector notation, our system can be expressed as 


N 

- cf = 0 , 

i=o 

for e = 1, 2,..., iVg and i = 0,1,..., iV. The local element {N + 1) x {N + 1) Galerkin 
differentiation matrix, is defined as 

Dlf = z,j=0,l,...,N. (2.68) 

Using LGL quadrature, can be calculated with the relationship 

Q 

= = = i,j = 0, (2.69) 

A:=0 


where the LGL weights given by Equation (2.49) and A is the Legendre PS 

differentiation matrix (2.57). Since LGL quadrature rule is exact for polynomial integrands 
of degree less than or equal to 2N — 1, the numerical integration is done exactly when 
Q = N LGL integration points are used. If Q = iV LGL quadrature nodes are used, the 
approximation to the {N + 1) x 1 RHS vector simplifies to 



-(e)m 




f{x 




for e = 1,2,..., iVe and i = 0,1,..., iV. 

Remark 2.13. So far, the required objects have been identified to solve the system numer¬ 
ically with element-based Galerkin. However, since nodal basis functions are continuous 
across element boundaries and LGL nodes include both endpoints, a global solution to our 
problem can be found. To do this, a global assembly or direct stiffness summation can be 
done, where the direct stiffness summation operator is [44] 
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The global equations to the problem become 

Np 

= 0, 1 = 1,..., Np. (2.70) 

j=i 

The global Np x Np Galerkin differentiation matrix, Djj, and RHS vector, , are then 
defined by 

Ne Ne 

Djj = /\Dlf and c^p^ = 

e=l e=l 

where Np = {NeN + 1) is the total number of grid points. Note that the direct stiffness 
summation operator, Afii’the mapping (i, e), (j, e) —)■ /, J [44]. So for the local 
differentiation matrix and RHS vector 

“oo “oi %N 

Uio Ujj 

U-ATO “ATI ^NN 

the direct stiffness summation operations Djj = f\ D^f and c^p^ = A result in the 

e=l e=l 

global Np X Np differentiation matrix and iVp x 1 RHS vector. 
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In order to guarantee feasibility of the discretized problem, the following inequality 


constraint is suggested, 


l\p 

- o'"-' 

J=0 


< 6^^, I 


0,1,...,Np, 


where 5^^ is the feasibility tolerance that is dependent on Np and the smoothness of x and 
u (see Section 5.2). The initial conditions and endpoint conditions may be approximated 
similarly by 


- Xo\\^ < and ||^ < 

Note that although the discretization for Problem (2.67) is element-based, the con¬ 
tinuous Galerkin formulation (2.70) is global in nature. This, however, is not the case for 
the discontinuous element-based Galerkin approach. 

23.2.2. Discontinuous Galerkin 

Consider a discontinuous element-based Galerkin approach to discretizing (2.52). 
Again, for this discussion, the LGL node structure will be used for the problem discretiza¬ 
tion. In this approximation, the weak integral form of (2.52) inside each element, takes 
the form [44] 


[ (i;(")^(t) - /(a;(")^(t), u(")^(t))) dt = 0, 

J Q,, 

for e = 1, 2, ..., Ae and i = 0,1, ..., A, where Q = Ijf^i defines the total domain. 
The state trajectory, x{t), is approximated inside each element, by interpolating A-th 
order Lagrange polynomials, at the nodes by the relationship 

N 

j=0 
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for e = 1, 2,, iVe, where the LGL nodes, mapped back to the 

physical space inside each element, fie- Also, let u^{t) be an interpolating function of 


N 


U 






i=o 


where of continuous functions (not necessarily polynomials) with 


the property = <5^. Therefore ^{e)Nj ^ for e = 1,2,..., A'e and 


(e)^ 


j = 0,1,..., A^, and similarly, The relationship between the physical 


time domain, t G [to, t/] = 
by [44] 


Al) ANe) 
'-0 ) ''AT 


, and the computational space, { e [—1,1], is given 




1 and = 


At('=) 


dt, 


and conversely, 

t = (^ + 1) + and dt = -^d^, 

where At*^®^ = is the size of each element, fie, which can be nonuniform in length. 

The Lagrange polynomial defined on the LGL computational domain is given by 

j=oyA ?i) 

j¥=i 

The state trajectory, x, can now be approximated inside each element, fie, by 


N 


j=0 
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where {0 ^(.^)}^q are the Lagrange polynomials defined on the LGL grid. Likewise, u^{^) 
is given by 

N 

j=0 


where = ^ij- In the computational domain, the system becomes 






'-1 


for e = 1, 2,..., iVg and i = 0,1,..., iV. Integration by parts on the first term yields the 
weak form relationship 




1 Ay-(®) 




for e = 1, 2,..., A"e and i = 0,1,..., A^, and in terms of our approximating polynomials 
we have 



N 




1 

-1 



2 





Remark 2.14. With the discontinuous element-based Galerkin approach, we let x, u and 
the basis functions be discontinuous across element edges. A numerical flux term acts 
as a jump condition between elements [44], Here, we consider the centered flux relation¬ 
ship, = I proposed by Delfour et al. [66], where e and q denote the 

element and its neighbor, respectively. 


50 



Integrating by parts, yet again, results in the Galerkin strong form relationship 







'-1 


for e = 1, 2,..., iVe and i = 0,1,..., iV. Since LGL nodes are used, the boundary term, 
7]^^\ may be simplified as 



I i , z = N, 

[o, 


(Ne) 

Vi = 


1 ^^(Ne)m _ ^iNe-l)NN^ ^ i = 0, 

0 , * 7 ^ 0 , 


for elements f2e = and respectively, and for each other element (f2e 7^ 
we have 


(e) / 

Vi = <, 


1 ^^(e)Aro _ ^ie-l)NN'^ ^ 
1 ^^(e+l)Af0 _ ^ie)NN'j 

0 , 


i = 0, 
i = N, 
i^0,N. 


In matrix-vector notation, our system may be expressed as 


N 

E 

j=0 


D 


(e) -(e) 


O 3 


xy + Z], 


(e) 


- 4'’ = 0, 


for e = 1,2,...,iVe and i = 0,1,... ,N, where the local element {N + 1) x {N + 1) 
Galerkin differentiation matrix, is the same as that defined in (2.68) and (2.69). If 
Q = N LGL quadrature nodes are used, the approximation to the (iV -(-1) x 1 RHS vector 
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simplifies to 



-(e)Ni 






for e = 1, 2,..., iVe and i = 0,1,..., iV. The loeal diseontinuous formulation therefore 
beeomes 

N 

= 0, (2.71) 

j=0 

for e = 1,2,..., iVe and i = 0,1,..., iV. In order to guarantee feasibility of the diseretized 
problem, the following inequality eonstraint is suggested, 


N 

E 

J=0 






-{e)Ni 


< 5 


N 


for e = 1, 2,..., iVe and i = 0,1,..., iV, where 5^ is the feasibility toleranee that is 
dependent on N and the smoothness of x and u (see Seetion 5.3). The initial eonditions 
and endpoint eonditions may be approximated similarly by 


^ OQ — ^ ' OQ — 


Note that unlike the eontinuous element-based Galerkin approaeh that ean easily be 
formulated for a global solution (2.70), the diseontinuous Galerkin formulation (2.71) is 
purely loeal in nature. The eommunieation between elements is done only by the boundary 
term, It is therefore easy to see that the diseontinuous Galerkin formulation is easy to 
parallelize for eomputational effieieney. Additionally, the flexibility and diseontinuous na¬ 
ture of the formulation lends itself to problems with eomplex geometries and diseontinuous 
solutions. 

Global speetral method teehniques, speeifieally, eolloeation (or PS methods) have 
beeome the method of ehoiee for diseretizing system dynamies in a number of direet meth- 
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ods for optimal control, such as the Legendre PS method. Element-based collocation tech¬ 
niques have also been investigated for the use in direct methods for optimal control by Ross 
et al. [12]. These methods have gained attention due to their flexibility as well as compu¬ 
tational efficiency. In Chapter 5, we will further investigate the use of the element-based 
Galerkin formulations for optimal control. However, we will first consider additional mo¬ 
tivation for the use of the weak integral formulation in Chapter 3. This will lead to the 
creation of a family of Galerkin-based formulations called, Galerkin optimal control. 
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CHAPTER 3: 

MOTIVATION FOR GALERKIN OPTIMAL CONTROL 


PS methods for optimal control have been shown to be good all-round methods for 
solving nonlinear optimal control problems. In particular, the Legendre PS method has 
gained much attention in recent years. As mentioned previously, two highlights are the 
successful use of Legendre PS method for the first ever zero-propellant attitude maneuver 
of the International Space Station [4] and the first ever minimum-time rotational maneuver 
performed in orbit by a NASA space telescope called TRACE [8]. The formulation of the 
Legendre PS method is provided in Section 3.2, but first consider the following problem of 
optimal control. 

3.1. Problem B (Bolza Problem) 

Determine the state-control function pair, t h-)■ {x, u) G x that minimizes 
the cost functional 

~ J F{x{t),u{t))dt + E{x{—l),x{l)), (3.1) 

subject to the dynamics 




(3.2) 


endpoint conditions 


e{x{-l),x{l)) = 0, 


(3.3) 
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and mixed state-control path conditions 


h{x(t),u(t)) < 0. (3.4) 

It is assumed that F : x R, E : R^^ x R^^ — )■ M, / : R^^ x — )■ R^^, 

e : R^^ X R^^ —)■ h : R^^ x —)■ are Lipschitz continuous with respect to 
their argument. It is also assumed that an optimal solution (a;*(-), «*(■)) exists. Additional 
assumptions related to the smoothness of a;*(-) and u*{-) are provided in the feasibility and 
consistency theorems included in Chapters 4, 5 and 6. 

3.2. Legendre Pseudospectral Method 

In the Legendre PS method approximation to Problem B, the states are approxi¬ 
mated with globally interpolating A^-th order Lagrange polynomials defined on LGL grid, 
{tj}f=o- Recall that the LGL points are defined by to = —1 < fi < ■ ■ ■ < = 1 and are 

the roots of Equation (2.42). The state trajectory, x{t), is approximated by 

N 

x(t) ^ it) = (t)x^^. 


x^^^x{tj), j = 0, 

and similarly, ^ u{tj). The Lagrange polynomials, of order N, are given by 

Equation (2.32), and have the property, = 6ij, for i, j = 0,1,..., At. 

In the Eegendre PS method, a solution to the differential equation x — f{x, u) = 0 
may be approximated at the EGE nodes with the following formulation 

N 

* = 0,1,..., At, (3.5) 

j=0 
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since, the derivative of x^{t) at each LGL point exactly equal to 

N 

x^iu) = 

i=o 

for any polynomial with degree less than or equal to N [46], where A is the Legendre 
PS differentiation matrix (2.57). A feasible solution to the equality dynamical constraint 
may not exist. In order to guarantee feasibility of the discretized problem, Gong et al. [5], 
suggest the following inequality constraint, 

N 

- f{x^\ <S^, t = 0,l,...,N. 

oo 

Remark 3.1. Note that 6^ is the feasibility tolerance and is dependent on N and the 
smoothness ofx and u. For x G (see Appendix C), m > 2 and u G C"^[—1,1], it has 

been proven by Gong et al. [5] that 5^ = (N — 1) 

The endpoint conditions and path constraints are approximated similarly by 

h{x^\u^^) <5^-1, i = 0,l,...,N, 

where 1 denotes [1,..., 1]^. Lastly, the cost functional J[x(-),u(-)] is approximated by 
LGL quadrature rule, 

N 

J[x{-),u{-)] ^ J^{x^,u^) = + 

i=0 

where x^ = x^^, ..., x^^] , ..., and the LGL 

weights (2.49) associated with the LGL points, To allow for a practical search area 
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for the optimal solution the following constraints are added 

G X, e [/, i = 0,1,..., N}, 

where X and U are the search regions that contain the optimal solution of the discretized 
nonlinear optimization. 

The resulting optimization problem can be solved using existing NLP algorithms. A 
feasible solution can be found that satisfies the tolerances specified in the NLP by adjusting 
the order of polynomial used in the approximation. The theoretical underpinnings of the 
Legendre PS method have been studied in great detail over the last two decades. Theorems 
for feasibility, consistency and convergence of the Legendre PS method approximations 
can be found in [3, 5, 6, 67, 68]. Although the Legendre PS method has been shown to 
produce accurate solutions on a wide variety of optimal control problems, it has proven to 
be a challenging task to modify this method to efficiently solve multi-scale problems, one 
for which the state(s) and control(s) evolve at different timescales. An example of such a 
problem is given next. 

Example 3.1. Consider the following boundary value problem given by Williams [69] of 
minimizing the cost function 



subject to the dynamics 

xi(t) = X 2 and X 2 (t) = C sm{kt) + u, (3.7) 

and with boundary conditions 

a;i(0) = 0, 0 : 2 ( 0 ) = 0, xi{tf) = 1 and X 2 {tf) = 0. (3.8) 


58 



The analytic solution to this problem is given by 


,, c . e c 

xi[t) = sin(/ct) - Cl — - C2— + —t, 

b 2 k 

,, c ,, , e c 

X 2 [t) = cos(fct) - Cl — - Cat + —, 
k 2 k 


u(t) = -Cit - C2, 


(3.9) 

(3.10) 

(3.11) 


obtained via Pontryagin’s maximum principle. The constants are defined as 

C fC C\ \2 f C C \ 

Cl = i^-cos{ktf) - ^ - -^tfj , (3.12) 

C2 = ^ (^^cosiktf) - (^1 + , (3.13) 

tf = 10, C = 0.1 and k = 8. This problem was solved using the Legendre PS method 
with optimality and feasibility tolerances of 5 x 10“^ and 5 x 10“^, respectively. The 
exact solution was used as an initial guess. Figure 12 shows the Legendre PS method 
approximations of order, iV = 50. 
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Figure 12: Exact solution and Legendre PS method approximation with N = 50 for Exam¬ 
ple 3.1. 


From Figure 12, it is apparent that xi and X 2 evolve on different timescales. This 
is confirmed by viewing the Legendre spectral coefficients of x\ and X 2 presented in Fig¬ 
ure 13. Note the difference in magnitude of the Xi and X 2 Legendre spectral coefficients, 
particularly between n = 5 and 40. 
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(a) Legendre spectral coefficients for state xi. 
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(c) Legendre spectral coefficients for control u. 


Figure 13: Legendre speetral eoefficients for xi, X 2 and u for Example 3.1. 


The differenee in evolution of the xi and X 2 system dynamies suggests that problem 
(3.6)-(3.8) may be approximated more effieiently using a multi-seale numerieal teehnique, 
where slow state, xi, and fast state, X 2 are diseretized on different timeseales. Consider the 
following general multi-scale optimal control problem, in which the slow and fast states. 
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Xs(t) and X fit), are associated with the slow and fast dynamics, respectively. This modified 
Problem B is presented as Problem B. 

3.3. Problem B (Multi-scale Bolza Problem) 

Problem B. Determine the state-control function, 1 1 — )■ ixs,Xf,u) G x x 
that minimizes the cost functional 

J[xsi-),Xfi-),ui-)] = j Fixsit),Xfit),uit))dt + EiXsi-l),Xsil),Xf{-l),Xfil)), 

subject to the dynamics. 


Xsit) = fiXsit),Xfit),uit)), 
Xfit) = gixsit),Xf{t),u{t)), 


endpoint conditions. 


e(a;,(-l),a;,(l),a;/(-l),a;/(l)) = 0, 

and mixed state-control path conditions. 


hixsit),Xfit),uit)) < 0. 


It is assumed that F: x x —)■ M, E: 


iN.s X X M, e: 




iN.s X X X ^ M, /: 


X X X 




,h: 


X R^^f X -)■ 


R^h 

are Lipschitz continuous with respect to their argument. It is also assumed that an 
optimal solution iKi-):X}i-),u*i-)) exists. 


A number of methods have been investigated for solving multi-scale problems 
such as Problem B, by casting the slow and fast dynamics of the problem onto differ¬ 
ent timescales. Recently, Desai et al. [70] and Williams [69] provided varied techniques. 
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While Desai et al. use an elemental approach where fast and slow dynamics are treated 
with similar order polynomials within different size subintervals and Williams uses a tech¬ 
nique where the slow dynamics are approximated with a weak formulation. Additionally, 
in [71], Gong et al. investigate the use of a Tau-like method to discretize the slow dynamics, 
after discounting a straightforward modified Legendre PS method approach. This modified 
Legendre PS method will be presented next for discussion purposes. 

3.4. A Modified Legendre PS Method for Multi-scale Problems 

Consider the following modified Legendre PS method approach to solving Prob¬ 
lem B. The states and controls are approximated with globally interpolating Lagrange 
polynomials on different LGL timescales. The slow state, Xs{t), is approximated on sparse 
grid while the fast state, Xf{t), on dense grid where M < N. The slow 

and fast states are defined by the following approximating polynomials 

M 

Xsit) ^ xfit) = '^(t)f{t)xf\ 
j=0 
N 

Xf{t) ^ X^{t) = 

j=0 

where the Lagrange polynomials {0^(t)}^o defined on grids {rj}^o 

and respectively. Let 

xJ^^Xfitj), j = 0,1,...,A^, 

and similarly, ^ u{tj), for j = 0,1,... ,N. 

Remark 3.2. For simplicity, the control variable, u(t), is approximated on the dense grid 
however this need not be the case. Modifications may be made to this method to 
cast the control onto a unique grid, such as sparse grid where M < N [71]. 
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A solution to the differential equations Xg = f{xs,Xf,u) and = g{xs,Xf,u) may 
be approximated by discretizing the slow dynamics over the dense grid with the following 
formulation 


M 

i=o 

N 

E aNN-NJ 
Xj 

3=0 


f{x^\xf,u^^) = Q, I 
g{x^\xf,u^^) = Q, I 


0 , 1 ,. 

0,1,...,iV, 


where is the standard (A^ + 1) x (A^ + 1) Legendre PS differentiation matrix (2.57) 
and is the (A^ + 1) x (M + 1) Legendre PS differentiation transformation matrix 
(2.39). The slow state approximation projected to the dense grid, , may be calculated by 
the linear mapping with the relationship 


j=0 


for i = 0,1,..., A^, where is the (A^ + 1) x (M + 1) transformation matrix (2.37). 

Remark 3.3. Projecting the slow dynamics onto the dense grid provides a way of capturing 
the high frequency information of the fast state [71]. 

The dynamical constraints therefore become 

II M II 




3=0 


<5^, i = 0,l,...,Ar, 


N 

3=0 


-g{xT.xf,u^^) 


< 6 ^, i = 0,l,...,Ar, 


64 



where 5^ is the feasibility tolerance. The endpoint conditions and path constraints are 
approximated similarly by 


\e{xr,xr,xr,xr)\L<s\ 

h{x^\xf\u^^) <6^-1, i = 0,l,. 


.,N. 


Lastly, the cost functional J[x{-),u{-)] is approximated by the LGL quadrature rule, 


N 


J[a;(-),M(-)] ~ = ^F{xf\x^\u^^)wi +E{xf^,xf^,xf^,xf^), 


i=0 


where the LGL weights (2.49) associated with the LGL points, {L}i=o 




Xf , Xj , 


r-MO -Ml -N ^ r 

[■^s , Jjg ■} ■ ■ ■ : J j •i'f [• 

and = [u^^, ..., 


- r^A'O -7V1 -NNl 


.,Xf j 


To allow for a practical search area for the optimal solution the following constraints are 
included: xf exs,xf e X f and e U, where Xg, X f and U are the search regions 
that contain the optimal solution of the discretized nonlinear optimization. 

Example 3.1 (continued). Consider again problem (3.6)-(3.8) solved with the proposed 
multi-scale Legendre PS method. The following analysis follows that given by Gong et 
al. in [71]. Here we discretize the slow state, xi, on LGL grid, {rj}^o, and fast state, 
X 2 and control, u, on LGL grid that This problem was solved 

with optimality and feasibility tolerances of 5 x 10“"^ & 5 x 10“^, respectively. The exact 
solution was used as an initial guess. Figure 14 shows the visual accuracy of the multi-scale 
Legendre PS method approximations with = 40, = 50 and iV„ = 50. A decrease 

in the approximation order of the slow state by 10 causes a significant decrease in the 
accuracy of the overall approximation. This is particularly apparent in the approximation 
of the control, u, in Figure 14. 
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Figure 14: Exact solution and multi-scale Legendre PS method approximation with = 
40, Nx 2 = 50 and Nu = 50 for Example 3.1. 


Also note that if now the control, u, is cast on a unique LGL grid such that Nu < 

Nx 2 , the NLP becomes infeasible. However, if u is cast on a unique LGL grid such that 
Nu > Nx 2 , an accurate solution is obtained. 

Although the Legendre PS method has been shown to produce accurate solutions on 
a wide variety of optimal control problems, approximating the derivative of a function using 
a standard PS differentiation matrix may introduce errors into the approximation. This may 
be an issue when using a multi-scale approach such as the one presented in Section 3.4. It 
should be mentioned, however, that the Tau-like method of Gong et al. [71] produces accu¬ 
rate solutions for this multi-sc ale approximation. However, to understand what happened 
with the straightforward multi-scale approach, we look to Jackson’s Theorem. 
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3.5. Jackson’s Theorem 


Jackson’s Theorem allows for a bounding of the spectral coefficients (2.26), 
of a function, H{t), in terms of the function’s derivative, h{t). 

Lemma 3.1 (Jackson’s Theorem). [72] Let h{^) be of bounded variation in [—1,1], Define 

= H{-1) + j\{s)ds, 

then the sequence of the spectral coefficients of H{f) satisfies the following in¬ 

equality 

|a„| < + 

for n > 1 where U{h{f)) is the least upper bound of |/i(OI V{h{f)) is the total 
variation of h{f) (see Appendix A). 

To see how this theorem affects a PS approximation of a function’s derivative, let 
H be the approximating polynomial error of a function, let h be its derivative, and let 
{an}'^=N+i be the spectral coefficients of H. Jackson’s Theorem says that even though 

OO 

\an\ may be very small (such as in the tail of the spectral coefficients dropped from 

n=N-\-l 

an approximation) the error in the approximation of the derivative may be relatively large. 


< {U{h{0) + V{hm. 

0 


This factor of could potentially add unnecessary errors when approximating a system’s 
dynamics using a standard differentiation matrix. 


Remark 3.4. This idea is further understood by considering the following estimates on the 
approximation of any function fit) G (see Appendix C). Consider if (f) approximated 
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by the truncated Legendre series 


N 

P^{t) = 

j=0 

where Lj is the Legendre polynomial of order j and {aj}^o spectral coefficients of(. 
The error estimate between ((t) and its approximation, p^(t), is given by 

oo oo oo 

iic(t)—p (^)ii^oo = ^ i%i’ 

j=N+l ^oo j=N+l j=N+l 

due to the property of the Legendre polynomials [46], \Lj{t)\ < 1, t G [—1,1] (see Ap¬ 
pendix A for definition of L^-norm). Additionally, the following estimate is provided 
by [46] 

IICW-P^WIL* SCiCoiV-i, (3.14) 

where Ci is a constant independent of N and Cq = Vthe total variation (see 

Appendix A). Now consider the error estimate between ((t) and its approximation 

N 

P^{t) = 

j=0 

given by 

OO OO ^ oo 

c(t)-p^(t) = r - 9 + 

L°° L°° Z 

j=N+\ J^oa j=N+l j=N+l 

due to the property of the Legendre polynomials [46], Lj(t) < lj(j + 1), t G [—1,1]. 
Additionally, the following estimate is provided by [46] 

m-p^(t) <CsC2Nf (3.15) 

L°° 
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where C 3 is a constant independent of N and C 2 = |C|/f 2 ;iv, the Sobolev seminorm of 
C (defined in Appendix C). The disparity between estimates (3.14) and (3.15) are clear 
and thus p^(t) and p^(t) may converge at different rates. In fact, while estimate (3.14) 
proves convergence of p^(f), estimate (3.15) shows that p^(f) may not converge at all. 
This disproportionate convergence behavior may have compounding effects when using 
the multi-scale approach for optimal control (introduced in Section 3.4). 

Example 3.1 (continued). Consider again problem (3.6)-(3.8) solved with the proposed 
multi-scale Legendre PS method and = 40, Nx 2 = 50 and Nu = 50. If now the 
optimality and feasibility tolerances are relaxed and decreased to 5 x 10-2 and 5 X io-\ 
respectively, the NLP constraints are satisfied and an accurate approximation of the states 
and control is obtained. In the context of Jackson’s Theorem and Remark 3.4, this should 
not be a surprise. From Figure 13, the Legendre spectrum of the dropped xi modes consist 
of coefficients with magnitudes of O (10“^). In fact with the lower optimality and feasibility 
tolerances, the multi-scale Legendre PS method can now produce accurate solutions for 
lower order approximations of xi, such as = 10 (with = 50). However, if 

a reduction in the control approximation is also the goal such as, < 40 (with = 10 
and = 50), a further reduction in the optimality and feasibility tolerances are required 
in order to satisfy the NLP constraints. 

The consequences of Jackson’s Theorem on multi-scale PS methods for optimal 
control are significant. Certainly, the class of problems that can obtain an advantage from 
this approach is limited. In general we can only hope to benefit from multi-scale PS when 
reducing the polynomial order of system variables that have extremely small Legendre 
expansion coefficients at the tail of the spectrum. This will require us to use a different 
approach if we hope to target a larger class of optimal control problems. Proposition 3.1 
highlights the advantage of an alternate method of discretizing the system dynamics, one 
in which the derivative of higher order terms does not disproportionally add to the overall 
error of the approximation. 
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3.6. Motivation for the Weak Integral Formulation 


Proposition 3.1. Let Lj(t) be the Legendre polynomial of order j. Suppose 


OQ 

j=N^l 

OQ 

= Y 

j=N+l 


are both uniformly convergent on [—1,1], then 


/I oo 

Lj(t)e'^(t)dt = 2 ttj, 

•1 AT I 1 


i-\-j odd 


(3.16) 

(3.17) 


(3.18) 


for all 0 < i < N. 
Proof 


/I „i 

Lj(t)e^(t)dt = Oj / Li(t)Lj{t)dt 

■1 j=N+l 

CO , „i 

= Y 1-1 - / Li{t)Lj{t)dt 

— Ar_Li V ./ — I 


j=Af+l 


Since the order of eaeh {Li(t)}^Q is less than N and the order of eaeh {7^j(t)}?^+i is bigger 
than N, the orthogonality of the Legendre polynomials implies 



Li(t)Lj(t)dt 


0 . 


Due to the following properties of the Legendre polynomial, 


Lfc(l) = 1 and Lfc(-l) = (-l)^ 


Equation (3.18) follows. 


□ 
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Remark 3.5. A similar result is found for the case that Lagrange interpolation polynomials 
are used as test functions. Let be the Lagrange polynomials of order N, defined 

on grid to = —1 < ti,... ,1^-1 < In = 1- Also, let Lj(t) be the Legendre polynomial of 
order j; letrjit) and pit) be defined by (3.16) and (3.17), respectively. Then 

/ I nl 

(f)r](t)dt = (j)f{t)Lj{t)dt 

■1 j=N+l 

= ^ Oj ((j)f{t)Lj{t) \-i- [ ^f{t)Lj{t)dt^ . 

j=N+l ^ ^ 

Since the order of the polynomials {0f^(t)}^o the order of each 

is bigger than N, the orthogonality of the Legendre polynomials implies 

j ^4>fit)Ljit)dt = 0. 

Due to the following properties of the Legendre and Lagrange polynomials, 

0f(4) = 4*, i>fc(l) = l and Lfc(-l) = (-l)^ 


we have 



{i)i"(t)dt = { 


j=A^+l 

oo 

j=N+l 


i = 0, 
i = N, 


0 , 


i^0,N. 
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Remark 3.6. The weak integral approximation to the derivative x(t) follows from Re 
mark 3.5. Let 


N 

j=0 

N 

j=0 


Multiplying x by a test function if) then integrating over the domain gives 



4>f {t)x{f)dt 



ff {f)x^ {f)dt + 



ff{t)e^{t)df 


for i = 0,1,..., N, where the residual term, e^(t), and its derivative, e^(t), can be ex 
pressed by 


OO 

j=N+l 

OO 

j=N+l 


From Remark 3.5, the weak integral residual term is 



( 


OO 


E %(-!)'+■. 

j=N+\ 


<t>f{t)Lj{t)dt = I 


OO 


j=N+l 


1 = 0 , 

i = N, 


0, i 0, N. 
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Therefore, the error in the weak integral differentiation term may be bounded as 


(j)f {t)x{t)dt 


'-1 


(t)dt 


Nt 


'-1 


( oo 

X/ II) 

j=N+l 


< i 


X/ II) 

j=N+l 

0 , 


i = 0, 

i = N, 


In other words, the accuracy of the weak integral approximation to x is related to the 
Legendre spectral coefficients of the dropped modes. If the Legendre spectral coefficients 
of the dropped modes are negligible, that is 


l«il < 

j=N+l 


where 5 1, then 



{t)x{t)dt 



{f)x^ {f)dt, 


for i = 0,1,..., iV. The weak integral approximation to x will be similar to the accuracy 
of the approximation to x itself since 



L°° 


j=N+l 


OO 

]foo j = N-\-l 


oo 

< \(^J 

j=N+l 


Jackson’s Theorem (Lemma 3.1) along with Remark 3.6 present persuasive argu¬ 
ments for the use of the weak integral approximation (a.k.a. Galerkin methods) in place of 
traditional collocation techniques for approximating system dynamics in direct methods for 
optimal control. We will see that Galerkin methods may be formulated to efficiently solve 
multi-scale problems (see Section 5.4). As a preview, Figure 15 shows a comparison of the 
exact solution to Example 3.1 with the multi-scale Galerkin optimal control formulation 
numerical solutions with = 3, = 43 and A„ = 1. This problem was solved with 
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optimality and feasibility tolerances of 5 x 10 ^ and 5 x 10 respectively, and the exact 
solution was used as an initial guess. 



Figure 15: Exact solution and GOCM-MS approximation with = 3, = 43 and 

Nu = 1 for Example 3.1. 

Simulations show that the multi-scale Eegendre PS formulation becomes infeasible for the 
multi-scale approximation with orders = 3, Nx 2 = 43 and = 1 for any reasonable 
set of optimality and feasibility tolerances selected. 

Although the multi-scale Galerkin optimal control formulation shows promise in 
solving multi-scale problems, the advantages of the weak integral form are not limited to 
this problem set. Additionally, Galerkin formulations allow for the weak imposition of 
boundary conditions. That is, end conditions may be enforced only up to the accuracy of 
the approximation itself. Remark 2.9 highlights this property. Galerkin formulations with 
weak enforcement of boundary conditions have been shown to produce improved accu¬ 
racies in many applications. A detailed discussion is given by Canute et al. (see Section 
3.7 of [46]). The Galerkin formulation with weak imposition of end conditions may also 
allow for problem discretizations with other than EGE points, such as EGR and EG (see 
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Sections 5.5.2.1 and 5.5.2.2, respectively). An important highlight of the Galerkin for¬ 
mulations is that the feasibility and consistency theorems are proved for problems with 
continuous and/or piecewise continuous controls (depending on the Galerkin formulation). 

Lastly, Galerkin methods, as shown in Section 2.3.2, may be easily formulated 
as element-based methods, both continuous and discontinuous (see Sections 5.2 and 5.3, 
respectively). These element-based formulations may have benefits in approximating so¬ 
lutions to optimal control problems with multiple stages or those with discontinuous solu¬ 
tions, such as bang-bang control problems. As compared to global methods, these element- 
based techniques may be formulated to require less computational effort and memory. Ad¬ 
ditionally, the discontinuous Galerkin formulation may advantage from parallel computing. 
Chapter 4 will introduce a new numerical technique for solving nonlinear optimal control 
problems founded upon the Galerkin methods outlined in Section 2.3. 
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CHAPTER 4: 

GENERAL GALERKIN OPTIMAL CONTROL EORMULATION 


There are four parts to the numerical solution to an optimal control problem using 
the direct method: discretization of the system dynamics, discretization of the state-control 
constraints, integration of the cost function and solving the NLP. In the Galerkin optimal 
control approach introduced in [73, 74], we use Galerkin techniques to discretize the sys¬ 
tem dynamics based on LGL quadrature nodes. Recall, that the LGL points, are 

the roots of Equation (2.42) and therefore include the endpoints, t = ±1. Thus the dis¬ 
cretization works in the interval of [—1,1] and will then provide the framework for our 
problem (e.g., the state-control constraints will be discretized at these nodes). Recall that 
LGL quadrature rule will provide zero error for polynomial integrands of less than or equal 
to 2N — 1 [45]. Linally, LGL quadrature rule will be used to integrate the cost function. 
The resulting optimization problem can be solved using existing NLP algorithms. 

In addition to the general Galerkin optimal control formulation, this chapter con¬ 
tains a number of important feasibility and consistency results. Theorems 4.1 and 4.2 prove 
that nonlinear program Problems GOCM-S and GOCM-S (presented in Section 4.2) have 
feasible solutions to Problem B, where controls may be piecewise continuous. Additionally, 
Theorems 4.3 and 4.4 prove that the general Galerkin optimal control numerical approx¬ 
imation is consistent. That is, nonlinear programming Problems GOCM-S and GOCM-S 
are consistent approximations to the continuous optimal control Problem B. 

4.1. Method for Approximation 

In the general Galerkin optimal control approximation to Problem B, the state tra¬ 
jectory, x{t), is approximated with globally interpolating A^-th order Lagrange polynomi- 
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als, {0^}^g, defined on a grid of LGL nodes, {tj}^^Q, 

N 

x{t) ^ x^{t) = . 

j=0 


Due to the property of the Lagrange polynomials, = Sij, we have 




X 


Also, let u^{t) be an interpolating function of 

N 

J-O 

where is any set of continuous functions (not necessarily polynomials) with the 

property = 6ij. In the general Galerkin optimal control approach, a solution to the 

differential equation x — f{x,u) = 0 may be approximated at the LGL nodes with the 
following weak integral formulation [44] 

/ (U^ ~ dt = 0, (4.1) 


for i = 0,1,..., iV. In terms of the approximating polynomials, the system of equations 
becomes 


N 


E 



dt 



<t>-f{x\ 


u^)dt = 0, 


for 1 = 0,1,..., iV, and in matrix-vector form is given by 

N 

'^DijX^^ -Ci = 0, i = 0,1,... ,N. 

j=0 
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The (TV + 1) X (iV + 1) Galerkin differentiation matrix, D, is defined by 

d(j)f (t) 

= j (4.2) 

and the (iV + 1) x 1 right-hand-side (RHS) vector, c, is defined as 

Ci = j 4>^{t)f{x^{t),u^{t))dt, i = 0,l,...,N, 

The Lagrange polynomials, }^ 0 ’ their derivatives, are given by defini¬ 

tions (2.32) and (2.38), respectively. If LGL quadrature rule is used with Q = N quadrature 
points, the differentiation matrix, D, can be calculated exactly by the relationship 

Dij = (4)-^(4)wfc = L J = 0,1,..., TV, (4.3) 

fc=0 

where the LGL weights, {wijfLQ, are defined by Equation (2.49) and A is the Legendre PS 
differentiation matrix (2.57). The RHS vector, c, may also be approximated with quadrature 
with the relationship 

Q 

Ci ~ '^4>^{tk)f{x^{tk),u^{tk))wk, (4.4) 

k=0 

for T = 0,1,..., iV. If again, LGL quadrature rule is used with Q = N quadrature points, 
the RHS vector approximation, , may be simplified as 

N 

c^" = ^(l>i{tk)f{x^{tk),u^{tk))wk = f{x^\u^")wi, i = 0,l,...,iV. (4.5) 

Remark 4.1. Recall that for LGL quadrature rule, integration is exact for polynomial 
integrands of degree less than or equal to 2N — 1. IfQ = (iV -f 1) integration points are 
used, the RHS vector will integrate exactly when f{x(t),u(t)) is linear in x(t) and u(t). In 
the case of a nonlinear function f, the accuracy of integration (and therefore the accuracy 
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of the overall approximation) can be improved by increasing the number of quadrature 
points Q. However, in most cases, increasing the accuracy of integration by increasing Q 
will significantly add to computation time due to the required interpolation of the state and 
control vectors. This will be discussed in greater detail in Section 5.5.1. 

Therefore system (4.1) may be simplified into the form 

N 

= 0, i = 0,l,...,N. (4.6) 

j=0 

Remark 4.2. Inform (4.6), the resulting Galerkin equations that must be satisfied are 

- f{x^\ j Wi = 0, (4.7) 

for i = 0,1,..., iV. This implies the following, 

N 

= (4.8) 

j=0 

for i = 0,1,..., iV. Note that (4.8) is the same set of equations that would be relaxed when 
using the Legendre PS method (Section 3.2). Hence, numerical solutions to system (4.8) 
found via the collocation method will satisfy the Galerkin relationships in (4.6). However, 
inequality versions of (4.7) and (4.8) are used for computational purposes. As suggested 
by Jackson’s Theorem, a solution of the inequality version of (4.7) may not satisfy (4.8). 
In fact, the analysis for the Galerkin numerical formulation is based on the L^-norm. As 
a result, the inequality bound for the Galerkin formulation is not simply a multiple of the 
quadrature weight. Shown in Equation (4.13), the upper bound of the inequality includes 
a factor of ^/w. However, this relationship draws a clear connection between the general 
Galerkin optimal control formulation and the Legendre PS method, and will be exploited 
in the proof of convergence (Theorem 4.4). 
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Remark 4.3. Due to the results of Gong etal. [5], we know a feasible solution to the equal¬ 
ity dynamical constraint may not exist. In order to guarantee feasibility of the discretized 
problem a relaxation of this constraint is used. 

The dynamical constraint becomes 

N 

i=o 

where 5^ is the feasibility tolerance. The endpoint conditions and path constraints are 
approximated similarly by 

h{x^\u^^) < 5^ ■ 1, i 

where 1 denotes [1,..., 1]^. Lastly, the cost functional J[x{-),u{-)] is approximated by the 
LGL quadrature rule, 

N 

J[x{-),u{-)] ^ J^{x^,u^) = + 

i=0 

where x^ = x^^, ..., x^^] and ..., . To allow for a practi¬ 

cal search area for the optimal solution the following constraints are added 

{x^^ e X, e [/, i = 0,1,..., N}, 

where X and U are the search regions that contain the optimal solution of the discretized 
nonlinear optimization. 


<5^, i = 0,l,...,N, 
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4.2. Computation Strategy 


The computation strategy for the Galerkin optimal control formulation with strong 
enforcement of BCs is presented in two forms. First, the strategy for the continuous prob¬ 
lem, in terms of the approximating polynomials is outlined, denoted as GOCM-S. Next, 
the discrete problem, discretized on a LGL grid is presented, denoted as GOCM-S. 

Definition 4.1. Function g{t) is called piecewise if 3 to = —1 < ti <■■■< tk = 1 
such that g{t) is on each subinterval (f*,fj+i), i = 0,... ,k — 1; limg{t), limg{t) and 
lim g{t) exist for i = 1,..., k — 1; and g{t) is either left or right continuous at each point 
U. 

4.2.1. Computation Strategy for GOCM-S 

The computational strategy of the GOCM-S is to find the feasible solution x^(t) e 
X and it) G U for the following cases: 

Case 1. u{-) is piecewise and x{-) G and piecewise C^, 

Case 2. m(-), i;(-) G (see Appendix C) and m>2, 

that minimizes 





F{x^ (t),u^ {t))dt + 1), a;'^(l)), 


(4.9) 


subject to the Galerkin constraints 




(t) {x^{t) - f{x^{t),u^{t))) dt 


< Mwf N 




i = 0,1,..., 


(4.10) 
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where a = ^ and (m — 1), for Case 1 and 2, respectively; M is a constant independent of 
N and 


{ h, h > 0, 

0 , h<0. 


(4.11) 


4.2.2. Computation Strategy for GOCM-S 

The computational strategy of the GOCM-S is to find the feasible solution eX 
and e U for the following cases: 

Case 1. u{-) is piecewise and x{-) G and piecewise C^, 

Case 2. m(-), x(-) e (see Appendix C),m > 2 and x^'^~^\t) is of bounded 

variation in t G [—1,1] (see Appendix A), 
that minimizes 

N 

J^{x^, u^) = 5^F(x^\ u^^)wi + E{x^^, x^^), (4.12) 

i=0 


subject to the Galerkin constraints 


N 


< Mwf N-^, i = 0,l,...,N, 

oo 

||e(x^°,x^^)||^ < MAT-", 

h{x^\u^^) < MN-^-1, i = 0,l,...,N, 


(4.13) 


where a = ^ and (m — 1), for Case 1 and 2, respectively; and M is a constant independent 
of A^. 


4.3. Feasibility of Solutions 

In order to guarantee feasibility of the discretization, Theorems 4.1 and 4.2 show 
that a relaxation of the dynamical equality constraint to inequality is required, for GOCM- 
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S and GOCM-S, respectively. However, first a buildup of lemmas are required for each 
theorem. 

Lemma 4.1. [46] Let p^(t) be the N-th order truncated Legendre series polynomial ap¬ 
proximation to C E if™, t G [—1,1], then 

||C(f) - IL2 < aiaoiV"™, V t e [-1,1] 

and 

||C(f) -p^(f)||ioo < aaOsiV^"™, Vt e [- 1 , 1 ], 

where ai and 03 are constants independent of N; oq = |Cliym;iv, the Sobolev seminorm of( 
(see Appendix C); 02 = the total variation (see Appendix A); and m > 0 . 

(p^(t)with the smallest norm ||C(f) — P^(f)l|^2 called the N-th order best polynomial 
approximation ofQ in the Lf-norm.) 

Lemma 4.2. Let ((t) = g(t) + hut^(t), t G [— 1 , 1 ], where (, Ut^ G H^, g G H^, and 
UtXl) = u(t — tc) is the unit step function defined by 

[0, -l<t<f„ 

= \ 

[1, t,<t<l. 

N 

Also, let p^(t) = 'YhVnLn be the N-th order truncated Legendre series polynomial ap- 

n=0 

proximation to C. Then 

||C(t) -p^(f)||i2 < biboN~^ + b2{to, h)N~f Wt G [- 1 , 1 ], L f- - 1,1 and \h\ < 00, 

where bi and 62 are constants independent of N, and bo = ||5'||//i. {p^(t) with the smallest 
norm ||C(f) — ||^2 called the N-th order best polynomial approximation of(( in the 

Lf-norm.) 
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Proof. Let 


N N 

9^ {i) = '^9nLn and = y^UnL^ 

n=0 n=0 

be the truncated Legendre series of g and ut^, respectively. Then 

N 

P^{t) = y^^PnLn = g^{t) + hUt^t) 

n=0 

N N 

~ ^ ^^ gnLn T h^^UnLn, 

n=0 n=0 

for t e [—1,1], where 

Pn = gn + hUn, U = 0, . . . , N. 

Therefore, 

N 

\\C{t) -P^{t)\\^2= C{t)-Y^PnLn 

n=0 1,2 

= 11 {9{t) - 9^{t)) + h (ut^t) - u^^{t)) 11^2 
<\\9{t) - 9^{t)\\L2 + 1^1 Iht,(t) -<(t) 11^2 

N N 

= 9{t) - y^,9nLn + \h\ Ut,{t) - y^UnLn 

n =0 l 2 n =0 l 2 

From Lemma 4.1, 

N 

gif) - y^,9nLn <hihQN~^, 

n =0 1,2 
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where is called the polynomial of best approximation of g in the L^-norm, bi is a 
constant independent of N and bo = ||5'||j:^i. Additionally, from [45], 


\h\ 


N 

n=0 


<\h\ 

L2 



(4.14) 


where the normalizing constants, {7n}^0’ the Legendre polynomials are given by 
Equation (2.27) and the Legendre expansion coefficients, {un}^=o^ are defined by Equa¬ 
tion (2.26). Due to the properties of the Legendre polynomials. 


Ln{t) = {L'n+lit) - L'n-iit)) , andL„(l) = 1, 


the Legendre coefficients have the relationship 

Wn = ^ / {Ln+l{t) - L'n-lit)) “ Ln+l{U)) , 

J tc 

or may be expressed by 


Wn\ — 2l(-^"--l(^c) ~ Ln+l{tc))\ < - {\Ln-l{tc) \ + | L„_l_l (tc) |) . 


Since the Legendre polynomial has the bound [72] 

4 (-)^ 

\Ln{t)\ < 1 1 , t 7 ^ — 1 , 1 , 

n2 (1 — 1^)4 

we have the following bound on Un, 



62 1 
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for n > 2 , where 


hitch) = 


\h\ 


712 (1 -t2)4 


r, 


This is since 


1 1 2 2 V 2 

-r “I-r ^-r — — n'>2. 

(n —1)2 (n + l)2 (n —1)2 


Therefore, Equation (4.14) has the bound, 


\h\ 


N 


Ut, 


(t) ^ ' ^UriLr 


n=0 


1 1 




L 2 


\n=N+l 


Kn=N+l 




However, from the Integral Test Estimate [75], 


Hence, 


1 ^ ^ /■“ 1 1 
lim / —dx <7 ^ < lim / —dx = 


\h\ 


N 


Ut. 


.(t) ^ ^fi nLr, 


i=0 


< 62-^ 2 , 


L 2 


and finally. 


({t) — p^{t) ^2 < biboN ^ + 62-^ S Vt G [—1,1], to 7^ —1,1 and \h\ < 00. 


□ 

Lemma 4.3 (Holder’s Inequality). [76] Let the Holder conjugates p and q be real numbers 
with the property that 


1 

P 


+ 


1 

q 


1, 
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where p > 1 and q > 1. Then for any arbitrary complex-valued sequences x = {^k}k=o 
and y = {vk}^^Q the following property holds 

N / N \p/N \q 

k=0 \k=0 / \k=0 / 

Moreover, when extended to integrals, Holder’s inequality takes the form 

I 1 

where f and g are assumed to be p-th and q-thpower summable, respectively, ont E [a, b]. 

Lemma 4.4. Let be the Lagrange interpolating polynomials of order, N, de¬ 
fined on LGL grid Then, there exists a positive integer, Nq, such that, for any 

N>No, 


| 0 ; 


N 

'i L2 


< pw^ < qN 2 ^ 


for each i = 0,1,..., iV, where {wi}^Q are the LGL quadrature weights associated with 
the LGL points, P q are positive constants independent of N. 

f N ^ 

Proof. From [46], the discrete norm, IICatHat = ( X] I'Cfcl 1 , has the property 

Vfc=o / 


\\^n\\l2 < P||'Ca''| AT) 


for ePn, where p is a positive constant, independent of N. Since G Pat, we have 


lkflL2<Hk* 


N\ 


I AT’ 
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for each i = 0,1,..., iV. Furthermore, from the property of the Lagrange polynomial, 
= Sij, we have 

1 

for i = 0,1,..., iV. Also, from [46] we have, for i = 1, 2,..., — 1, 



N 


(1 - (L)^)2 <VJi< 


-(1 


(L 


for constants, 0 < Ci < C 2 , independent of i and A^; for i = 0, N, we have 


Therefore, 


W0,W7V 


2 

N{N + 1)' 


for each i = 0,1,..., A^. Finally, 

||0f 11^2 < pwf < qN-'^ 

holds for all N > Nq, where g is a constant independent of i and N. □ 

4.3.1. Feasibility of GOCM-S 

Theorem 4.1 (Feasibility of GOCM-S). Given any feasible solution t h-)■ (x, u), for Prob¬ 
lem B, consider the following two cases: 

Case 1. u{-) is piecewise and x(-) G and piecewise C^, 

Case 2. u(-), x(-) G H'^~^ and m > 2. 

Then, there exists a positive integer Nq such that, for any N > Nq, GOCM-S has a poly- 
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nomial feasible solution, (t), (t)) such that 


\\x{t)-x^{t)\\^, < MN-^ and \\u{t) - {t)\\< MN-^, 

where a = ^ and {m — 1), for Case 1 and 2, respectively; and M is a positive constant 
independent of N. 

Proof Let p{t) be the (iV — l)-th order truncated Legendre polynomial approximation of 
x{t). By Lemmas 4.1 and 4.2 there is a constant cq independent of N, for any N > Nq, 
such that 


\x{t) -p(t)||i 2 < CqN 


where a = | and (m — 1), for Case 1 and 2, respectively. Define 


x^ (t) = J p{s)ds + x{—l). 


Thenp(t) = x^{t) and 


x{t) - x^{t )\\^2 < 2coiV ", 


since, from Holder’s inequality (Lemma 4.3), we have 


x{t) — x^{t) = 


{x{s) — p{s)) ds 


'-1 


< / |i;(s) — p(s)|(is 




<a/2 [ I |i;(s) - p(s)p(is ] = A/2||i;(t) - p(t)||^2 < v^coiV 


(4.15) 


Let u^{t) be the iV-th order Legendre polynomial so that 


u{t) - M"(t) ^2 < CiN ". 
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From our Galerkin approximation, Holder’s inequality (Lemma 4.3), and Lemma 4.4, we 
have for each A; = 0 , 1 ,..., iV, 

1 

j <Pk (t) dt 

-1 

1 

-1 

< C2W^\\x^{t) - 

< C2wl\\x{t) -X^{t)\\^2 +C2wl\\x{t) - f{x^{t),U^{t))\\^^ 

= C2wi\\x{t) - p{t)\\^2 + C2W^\\f{x{t),u{t)) - f{x^{t),U^{t))\\^^ 

= CoC2w|iV“" + C2liw^ ||a;(t) - a;^(t)||^2 + C 2 I 2 WI ||M(t) - M^(t)||^2 

< CPC 2 Wl iV“" + C0C2I1WI iV““ + C1C2I2WI iV“", 


where {wk}k=o LGL quadrature weights and h and I 2 are the Lipschitz constants of / 
with respect to x and u, respectively, which are independent of N. It follows that 


J if) {x^it) - fix^it),u^{t))) dt 


< MwlN-'^ 


1-1 


holds for each /c = 0,1,..., iV, and all N > Nq, where M is a constant independent of N. 
For the endpoint condition we have 


a;(l) — a;^(l) = 


(i;(s) — pis)) ds 


1-1 


< / |i;(s) — p(s)|(is 




<a/ 2 ( / |i;(s) - p(s)|^(is ) = A/2||i;(t) - p(t)||^2 < i/2coiV ", 
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so we have, by Lipsehitz eondition, 


e(a;'^(-l),a;'^(l)) < MiV 


For the path eonstraint let V = [t\h{x^{t),u^{t)) > O} ,r> = [—1,1] \V, sinee 
h{x{t),u{t)) < 0. Then 


\h^{x^{t),u^{t))\\^^ = ( / {h{x^{t),u^{t))fdt 

'v 


< ( I {h{x^ it)) — h{x{t),u{t))Ydt 


< 


{h{x^{t),u^{t)) — h{x{t),u{t))Ydt + / {h{x^{t),u^[t)) — h{x{t),u{t))Ydt 
{h{x^(t), (t)) — h{x{t),u{t))Ydt 


'V Jv 

'*1 \ 2 


'-1 


= \\h{x^(t),u^{t)) - h{x{t),u{t))\\^, 

< h\\x{t) -x^{t))\\^^ + u\\u{t) -M^(t))||^2 < MN~'^, 


where and U are the Lipsehitz eonstants of h with respeet to x and u, respeetively, whieh 
are independent of N. Henee 


Thus a solution {x’^{t), {t)) to GOCM-S is feasible! □ 

4.3.2. Feasibility of GOCM-S 

Theorem 4.2 (Feasibility of GOCM-S). Given any feasible solution t h-)■ {x, u), for Prob¬ 
lem B, consider the following two cases: 

Case 1. u{-) is piecewise and xf) G and piecewise C^, 

Case 2. «(•), x{-) G and m > 2. 

Then, there exists a positive integer Nq such that, for any N > Nq, GOCM-S has a feasible 
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solution, {x^,u^) such that 


where a = ^ and {m — 1), for Case 1 and 2, respectively; and M is a positive constant 
independent of N. Additionally, u^(ti) = u{ti),fori = 0,1,..., iV. 

Proof. Let p(t)be the {N — l)-th order truncated Legendre polynomial approximation of 
x{t) in the L^-norm. By Lemma 4.1 there is a constant di independent of N, for any 
N > Ni, such that 


WHt) -pWIIloo < diN 

where (3 = (^^ — |), for Case 2. For Case 1, we refer to [77-79] which show the truncated 
Legendre approximation for discontinuous functions with jump discontinuity (such as the 
step function defined in Lemma 4.2) displays Gibbs phenomenon. However, the maximum 
amplitude of the overshoot has a finite limit; we conclude that for Case 1, = 0. Also, by 

Lemma 4.2 there is a constant d 2 independent of N, for any N > N 2 , such that 

II^W -pWIIl2 < d2N-^, 

where a = ^ and (m — 1), for Case 1 and 2, respectively. Define 

x^{t) = j p{s)ds-\-x{—l). 

Thenp(t) = x^(f) and 


x{f) — x^{t) ^2 < 2 d 2 N ", 
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since, from Holder’s inequality (Lemma 4.3), we have 

\x{t) — {t)\ = j {x{s) — p{s)) ds < j |i;(s) — p(s)|(is (4.16) 

<V2 (^J |i;(s) - p(s)|^(is^ = V2\\x(t) - p(t )\\^2 < V 2 d 2 N~'^. (4.17) 

Also, let u^{t) be an interpolating function of u{t), 

N 

j=0 

where is any set of continuous functions (not necessarily polynomials) with the 

property = 6ij, and therefore 

= u{tj). 

Since x^{t) is a A^-th order polynomial, we have 

N 

x^{tk) = y^Akjx^\ 

3=0 

where A is the (A^ + 1) x (A^ + 1) Legendre PS differentiation matrix (2.57) and 

x^^ = x^(h). 

Recall that the LGL quadrature weights, {wkj^^Q, have the property 

Wk < d3N~\l - /c = 1,2,..., A^ - 1, 
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for constant (is > 0 independent of k and iV; for /c = 0 , iV, we have 


Wo,Wn = 


N{N + l)' 


So, we have 


N 




Nj _ -Nk 


j=0 
.N, 


N 




3=0 


Wk 


= (4) - f{x^^, u^'')\wk = |p(4) - f{x^'', \wk 

< b(4) - f{x{tk),u{tk))\wk + |/(a;(4),M(4) - f{x^^,u^^)\wk 
= b(4) - x{tk)\wk + |/(a;(4),M(4) - f{x^'",u^^)\wk 
< ||p(t) - x{t)\\^^Wk + /i|a;(4) - x^{tk)\wk 
< d4WkN~^ + V2lid2WkN~°', 


for each k = 0 , 1 , 2 ,..., iV, where /i is the Lipschitz constants of / with respect to 
Putting this all together, we conclude that 


N 




::Nk 


3=0 


< diWkN-^ + V2hd2WkN-'^ < MwlN-'^, 


for each A; = 0 ,l, 2 ,...,iV, and all N > N^, where M is a constant independent of N. 
For the endpoint condition, we have 


x{tN)—x^{tN) = a;(l) — (1) = / |i;(s) — p(s)|(is 


AT/ 


'-1 


r-1 

'-1 


<72 |i;(s) - p(s)|^(is ) = V2\\x{t) - p{t )\\^2 < V 2 d 2 N 


So, by Lipschitz condition. 


e{x^{to),x^(1^)) <MN 


X. 
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For the path constraint, the following estimate holds 


\h{x{t),u{t)) — < l 2 \x{t) — x^{t)\ < V 2 l 2 d 2 N ", 

for each k = 0,1, 2,..., iV, where I 2 is the Lipschitz constants of h with respect to x. 
Hence, 


h{x^\ < h{x{tk),u{tk)) + MiV-" ■ 1. 

Thus a solution u^) to GOCM-S is feasible! □ 

Remark 4.4. Although, Theorems 4.1 and 4.2 do not provide exact feasibility tolerances for 
the existence of solutions to the GOCM-S and GOCM-S, we can be confident that solutions 
do in fact exist. Precise bounds may be found experimentally, using a recursive refinement 
process, by increasing the order of the approximation, N, until all the constraints in the 
NLP are satisfied. 

4.4. Consistency of Solutions 

Theorems 4.1 and 4.2 show that solutions exist to the GOCM-S and GOCM-S, re¬ 
spectively. However, the question still remains—will these solutions converge to those that 
we seek? The answer is yes—Theorems 4.3 and 4.4 presented below show that solutions 
to GOCM-S and GOCM-S, will in fact converge to the optimal solution of Problem B. 
However, first a definition and lemma are required. 

4.4.1. Consistency of GOCM-S 

Definition 4.2. The orthogonal system is complete in L^, t e [—1,1], if and 

only if, for ,^(t) G the condition = 0, V/c > 0, implies ||'C(f)|lL 2 = 0. 
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Lemma 4.5. Let e^{t) G with m > ^ (see Appendix C). Assume 




Nf 


1-1 


< bow!5^, 


for each i = 0,1,..., N, where {0f^(t)}^o ^^-e Lagrange interpolation polynomial of 

order N defined on LGL grid, ^be associated LGL quadrature weights 

and ho is a constant independent of N. Also assume 3 e(t) G H'^ with m > \, so that 


e - e 


N\ 


Il2 


< hi5 


N 


where 5^ < N °‘,a > \ and hi is a constant independent of N. Then 


Ikllis - 0. 

Proof Recall from Equation (2.35) that the orthogonal Legendre polynomials, [Lj (t) }^ 0 ’ 
can be written as linear combinations of Lagrange polynomials, [fifif) }^q, defined on the 
LGL grid, {tk}k=o^ 

N N 

^=0 2=0 
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where V is the generalized Vandermonde matrix (2.30). From Holder’s inequality (Lemma 4.3) 
and Lemma 4.4, for each i = 0,1,..., iV, we have 


< 


4>f{t)e{t)dt 

c 

-1 




(t) (e^(t) - e^(t) + e{t)) dt 


(e(^) - dt 


+ 


(j)f {t)e^ {t)dt 


'-1 


< I Wit) (e(t) - e^(t)) \dt + bowfd^ 
'-1 


< ||0fW|L2|Kt) -€^{t)\\^^ + bowi5^ 


< b2wf 5^ + b{)W 2 5^'^ = 63 w/ (5 




From [46], we have 


1 = 1 , 2 , 


.N-1, 


for constants 0 < 64 < 65 , independent of i and iV; for i = 0, iV, we have 


Wo,Wn = 


N{N + 1)' 


Also recall, from [72], that 


\Ljit)\ < - 


4(: 


j2(l -t2)4 


r, t^-l,l. 
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Since 


f-i < b^wf5^, for each i 


0,1,..., iV, we have 


Lj{t)e{t)dt 




N 1 

i=0 


< 


TV 1 

/ T-f 

i =0 -^-1 


{t)e(t)dt 


0i it)e{t)dt 

TV 


< 




i=0 


< 5 


Af 


2^2 


AT-l 


= 5 


N 


(iV(iV + l))2 

2^2 




2 = 1 


1 


) 


ivi 

) 


iV - 1 

V+& 6 TT^) <M^M- 
(iV(iV+l))2 J2iV2 / Vj 


N\ 2 


for each j = 1,..., iV, and constant, bj, for all N > Nq. Since <5^ < iV " and a > |, it 
follows that when iV —)■ cx) we have 



Lj(t)e(t)dt 


0 , 


for each j = 0,1,..., iV. Since the Legendre polynomials, are complete in 

space [80] and e ± Lj for all j = 0,1,..., iV, we can conclude. 


e(t )||^2 = 0 . 


□ 

Theorem 4.3 (Consistency of GOCM-S). Suppose {t),u^{t)) is a solution ofGOCM- 
S and there exists {x(t),u(t)) such that u{-), x{-) G with m > 2. Also, suppose 

x^(t) —)■ x(t) uniformly, and 

Ht)-x^{t)\\^,<K5^, (4.18) 

Ht)-u^{t)\\^,<K6^, (4.19) 
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where 5^ < N a > ^ and K is a constant independent ofN. Then {x(t),u(t)) satisfies 

II^W - = 0, 

e(a;(-l),a;(l)) = 0, 
h{x{t),u{t)) < 0, 

and is an optimal solution to Problem B. 

Proof. Let e^{t) = x^(t) — f {x^(t), (t)) and e(t) = x{t) — f{x{t),u{t)). From 

Lemma 4.5, to prove || <t)IU. = o, it is enough to prove 

j <Pk <CoW^6^, 

for each A; = 0,1,..., iV, where < iV““, a > f . Consider 

< j (Pk (i) «^(^))) dt + j (t) {x{t) - x^{t)) dt 

+ j (t)k{i){f{x^{i)^u^{i))-f{x,u))dt 

< C2w 15^ + wl\\x{t) - x^+ wl\\f{x{t),u{t)) - f{x^ 

< C2w 15^ + wl\\x{f) - x^+ liwl\\x{t) - x^{f)\\^^ + l 2 wl\\u{f) - {f)\\^^ 

< C2w15^ + Kwl5^ + hKwld^ + hKwld^ < cowl5^, 

where li and I 2 are the Lipschitz constants of / with respect to x and u, respectively, which 
are independent of N. It follows, from Lemma 4.5, that as iV —)■ cx) we have 

l|e|lL 2 = \\x{t) - f{x{t),u{t ))\\^2 = 0. 
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For the endpoint condition, since x^{t) —)■ x{t) uniformly, we have 

—)■ a;(l) and a;^(— 1 ) —)■ a;(— 1 ). 

Since, from the formulation of the computational strategy we have 
we conclude that e(a;(—l),a;(l)) = 0 as iV —)■ cx). 

For the path constraint, since h{x{t),u{t)) is piecewise if h{x{t*),u{t*)) > 0, 
3 an interval (a, b) in which h{x{t),u{t)) > 0. Then 

||/i(a;(t),M(t))||^2(„;,) = >0. (4.20) 

However, 


\\Kx{t),u{t))\\^ 2 ^a,b) < \\Hxit):U{t)) - + \\h-^{x^{t),u^{t))\\ 

< \\h{x{t),u{t)) - h{x^{t),u^{t))\\^^^^^^^ + MN-^ 

< hWxit) - a;^(t))||^2 + U\\u{t) - M^(t))||^2 + 


L^{a,b) 


where (3 and U are the Lipschitz constants of h with respect to x and u, respectively, 
which are independent of N. Hence, this is a contradiction, therefore h{x{t), u{t)) < 0 as 
N ^ 00 . 

Suppose that {x{t),u{t)) is not optimal. Then 3 {x*{t),u*{t)), so that 


J (a;*(-), «*(■)) < J 


Also, 3 {x* (t), u* (t)) such that 


x*^{t) - x*{t)\\^, < MN-'^ and 
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where {x*^{t),u*^{t)) is a feasible trajectory of GOCM-S. Therefore 

J > J . (4.21) 

However, 


\J {x^{-),u^{-)) - J (a;(-),M(-))| 

< j \F{x^{t),u^{t)) - F{x{t),u{t))\dt + |E(a;^(-l),a;^(l)) - E(a;(-1),a;(l))| 
< V2\\F{x^{t),u^{t)) - F{x{t),u{t))dt\\^^ + |E(a;^(-l),a;^(l)) - E(a;(-1), a;(l)) 

Due to the Lipschitz condition and assumptions (4.18) and (4.19) we have 


lim J ix^i-) 

N^oo ^ ^ 




J {x{-),u{-)) 


0 . 


Similarly, 


lim 

N^oo 




0 . 


Therefore, from (4.21) we have 


J {x*{-),u*{-)) > J (a;(-), «(■))• 


This is a contradiction, since we assumed 


J {x*{-),u*{-)) < J (a;(-), «(■))• 


We conclude that {x{t),u{t)) achieves an optimal cost and therefore is an optimal solution 
to Problem B! □ 
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4.4.2. Consistency of GOCM-S 

Lemma 4.6 (Theorem 6.5.5, [81]). Let f be Riemann integrable in [—1,1] and {0^(t)}^Q 
be the Lagrange polynomials of order N defined on a LGL grid, Define 

N 

/"(«) = 

k=0 

where = f^{tk),for k = 0,1,..., N. Then 

^ pi pi 

lim ^f{tk)wk = lim / f^{f)dt = / f{f)dt, 

k=0 

where {wk}k=o quadrature weights associated with the LGL points, {tk}k=o- 

Theorem 4.4 (Consistency of GOCM-S). Suppose h^^), 0 <k< is a se¬ 
quence of solutions to GOCM-S, {t h-)■ {x^ {t),u^ are their interpolating func¬ 
tions and there exists functions {x(t),u(t)) such that uf), xf) G with m > 2 and 

j,(m-i)(^) of bounded variation in t G [—1,1] (see Appendix A). Also, suppose 

lim ||M(t) - u^{t)\\ = 0, (4.22) 

lim ||i;(t) - x^{t)\\ = 0. (4.23) 

Al i. rv-, ’' I > J-j 


Then {x{t)^u{t)) satisfies 

'' 

ll^(^) - = 0, 

< e(a;(-l),a;(l)) = 0, 
h{x{t),u{t)) < 0, 

and is an optimal solution to Problem B. 
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Proof. This proof follows the outline of Theorem 2 given by Gong et al. in [5]. First, from 
assumptions (4.22) and (4.23) it is easy to show 

lim \\x{t) - x^{t)\\ = 0. 

Next, suppose that {x{t),u{t)) is not a solution. Then there is a time r G [—1,1] such that 

^ 0 . 

From [82], the LGL nodes, dense when iV —)■ cx). Then there exists a sequence, 

where 0 < < iV, and with property 

lim tjN = r, 

N^oo 


SO that 


x{t) - f{x{T),u{T)) = lim (T^(tiiv) - /(a;^(fiiv),M^(fjiv))) ^ 0. (4.24) 

N^oo ' ' 

Also, we have 

N 

X^{tiN) = ^AiNjX^f 
j=0 

where A is the Legendre PS differentiation matrix (2.57). Therefore, from Theorem 4.2, 


N 


N 


= 


j=0 


i=o 


- f{x^{tiN),U^{tiN))\WiN = MWiNN^^ 
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where w^n, is the LGL quadrature weight associated with LGL point, t^N, for each . This 
implies 


lim {x^{tiN) — f (x^(tiN), (tiN))) = lim MN^ ™ = 0, 

which contradicts Equation (4.24). We conclude that {x{t),u{t)) is a solution. 

For the path constraint, we consider the same contradiction argument given above. 
For the endpoint condition, since x^{t) —)■ x{t) uniformly, we have 

—)■ a;(l) and a;^(— 1) —)■ a;(— 1). 

Since, from Theorem 4.2, we have 


we conclude that 


e(a;(—l),a;(l)) = lim e(a;^(—1),a;^(l)) = lim = 0. 

N^oo N^oo 


For the cost functional we have 


N 


NO ;^NN\ 


k=0 


and 


J{x{-),u{-)) = / F{x{t),u{t))dt + E{x{-l),x{l)). 


'-1 
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By Lemma 4.6, we have 


N 


'-1 


F{x(t),u(t))dt = lim y^F{x(tk),u(tk))wk, 

N^oo^-^ 


k=0 


and therefore 


F{x{t),u{t))dt 


'-1 


N 


N 


-Nk -Nk 


[ ) F{x^’^,u^^)wk + 22 [Fix{tk),u(tk) - F{x^^, 

k=0 k=0 


However, 


N 

lim V] [F{x{tk),u{tk)) - F{x^^,u^'^)] Wk 

k=0 
N 

< lim '^\F{x{tk),u{tk)) - F{x^{tk),u^{tk))\wk 

k=0 

N N 

< li lim y^\x(tk) - x^(tk)\wk + k lim y^\u(tk) - u^(tk)\wk = 0, 

A /— f ^ ' ’ /V — t ^ ' ’ 


k=0 


k=0 


where /i and I 2 are the Lipschitz constants of F with respect to x and u. Thus we conclude. 


i: 


N 


F{x{t),u{t))dt = lim 'S^F{x^’^,u^’^)wk. 


k=0 


Finally, by Lipschitz condition. 


lim E{x ^,x ) = 1),a;(l)), 

N^oo 


and the limit 


lim J^{x^,u^) = J {x{-),u{-)) 

N^oo 


(4.25) 


follows. 


Wk 
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Let {0 < k < be an optimal sequence of solutions to GOCM- 

S and {t !-)■ be their interpolating functions. From Theorem 4.2, 

u*^{t)\\^^=0, (4.26) 

x*^{t)\\^^=0, (4.27) 


lim \\u*(t) — 

7V-5>oo" 

lim ||a;*(t) — 
Ar^>oo" 


and from (4.25) we have 

J{x*{-U*{-)) < J{x{-)M-)) = lim < hm {x*^{-U*^{■)). 

N—>-oo N—>-oo 

Finally, from conditions (4.26) and (4.27) we conclude 

J (a;*(-), «*(■)) = J {x{-),u{-)). 

Hence {x{t),u{t)) achieves an optimal cost and therefore is an optimal solution to Prob¬ 
lem B! □ 

Remark 4.5. Theorems 4.3 and 4.4 provide confidence that solutions not only existence to 
the GOCM-S and GOCM-S, but solutions will converge to the optimal solution. However, 
questions still remain about the conditions under which assumptions (4.18), (4.19), (4.22) 
and (4.23) exist (as pointed out by Gong et al. [5]). Answers for similar questions have been 
provided for the Legendre PS method by Kang [6], but like analysis for Galerkin optimal 
control is above the scope of this dissertation. 
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CHAPTER 5: 

ALTERNATIVE EORMS OE GALERKIN OPTIMAL CONTROL 


An advantage in using Galerkin optimal control is that it is a versatile family of 
formulations. There are a number of Galerkin forms that can be used to suit the problem at 
hand. In Chapter 4, the general formulation, Galerkin optimal control with strong enforce¬ 
ment of end conditions, was presented. This serves as the first of three global formulations 
that are outlined in this dissertation. The second global formulation is Galerkin optimal 
control with weak enforcement of boundary conditions, and will be discussed in Section 
5.1. (Additionally, a third global Galerkin optimal control formulation with Legendre test 
functions will be presented in Chapter 6.) Important results in this chapter include Theo¬ 
rems 5.1 and 5.2, which prove that nonlinear program Problems GOCM-W and GOCM-W 
(outlined in Section 5.1.2) have feasible solutions to Problem B, where the controls may be 
piecewise continuous. 

Next, the element based formulations will be presented and are divided into two 
forms: Galerkin optimal control with element-based continuous and discontinuous Galerkin 
techniques, which will be presented in Sections 5.2 and 5.3, respectively. As alluded to in 
Chapter 3, the Galerkin weak integral form improves feasibility of the multi-scale approach 
highlighted. The method of approximation for the multi-scale Galerkin optimal control for¬ 
mulation will be outlined in Section 5.4. Finally, Section 5.5 will discuss modifications to 
the Galerkin optimal control formulations, such as over-integration of the RHS vector and 
the use of quadrature points other than LGL, such as LG and LGR points. 

5.1. Galerkin optimal control with Weak Boundary Condition Enforce¬ 
ment 

The general Galerkin optimal control strategy presented in Section 4.2 describes a 
formulation in which boundary conditions are enforced in a strong sense, via a constraint 
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of the form 


Recall, this boundary enforcement method is also incorporated into the Legendre PS method 
(see Section 3.2). This section presents an alternative formulation of Galerkin optimal con¬ 
trol introduced in [73, 83], one that allows for enforcement of the problem end conditions 
in a weak sense through the dynamical constraint. 

5.1.1. Method for Approximation 

We now consider the Galerkin optimal control formulation with weak enforcement 
of boundary conditions. In this approximation to Problem B, the state trajectory, x(t), 
is approximated with globally interpolating iV-th order Lagrange polynomials, 
defined on a grid of LGL nodes, 

N 

x{t) ^ x^{t) = 

j=0 

Due to the property of the Lagrange polynomials, 4>^{U) = Sij, we have 

x^^=x^{tj), j = 0,l,...,N. 

Also, let u^{t) be an interpolating function of 

N 

«"(i) = (ijfl"'. 

where is any set of continuous functions (not necessarily polynomials) with the 

property'^j^(L) = 5ij. As done previously, taking the weak integral form of i;—/(a;,M) = 0 
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yields [44] 


) dt = 0, 


'-1 


dt 


for i = 0,1,..., iV. Integration by parts on the first term results in Galerkin weak form, 


-■i 


d(f). 


dt 


* x^dt + [4>^x^Y_i ~ J f{x^■,u^)dt = 0 . 


In terms of our approximating polynomials and the true boundary eonditions (letting a;^(—1) —)■ 
a;(—1) and a;^(l) —)■ a;(l)) we have 





f{x ,u )dt = 0, 




for f = 0,1,..., iV. Integration by parts, yet again, results in Galerkin strong form with 
weak enforeement of BCs, 


N „i 
j=0 "'-1 

/ N N 


' ^ \ 

- a;(-l) J 

.i=o / 

-J (jyff{x^,u^)dt = 0. 


The expression may be simplified as 

N 

'^^DijX^^ + Ki — c^* = 0, 

j=0 


for eaeh f = where the Galerkin differentiation matrix, D, and the RHS veetor 

approximation, c, are unchanged from those given in the GOCM-S methodology—given 


111 



by Equations (4.3) and (4.5), respectively—and 


xNO _ a;( — 1), 

i = 0, 

a;(l) — x^^, 

i = N, 

0, 



The BC term k now provides a natural way to introduce end conditions into our numer¬ 
ical scheme. BCs such as e(a;(—1), a;(l)) = \x{—l)—x^,x{l)—x^Y' = [0,0]^ can be 
imposed by defining k as 


_ ^0 

i = 0, 


i = N, 

0, 



for f = 0,1,..., iV. 


Remark 5.1. A similar technique is discussed by Ross et al. in [11, 84, 85]. The ability 
to weakly enforce boundary conditions is not limited to the Galerkin formulation. The 
collocation method (in the context of the Legendre PS method) may be formulated for weak 
enforcement of end conditions by modifying the discrete differential Equation 3.5 with a 
boundary condition term k 


N 

+ ki- f{x^\u^") =0, f = 0,1,...,iV, 

j=0 

where k is given by 


ki 


Wi ’ 

^ xf-x^N 


(5.1) 


0 , 
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i = 0, 
i = N, 
i 



for fixed end conditions e{x{—l),x{l)) = [a;(—1) — — x^~\'^ = [0,0]^. 

Remark 5.2. For existing direct methods for optimal control it is common to enforce the 
problem’s BCs in a strong sense (or exactly) by making them a set of constraints enforced 
by the nonlinear program (NLP) [5]. With the Galerkin optimal control formulation with 
weak boundary enforcement, BCs can now be enforced in a weak sense. In other words, 
BCs can be imposed only up to the order of accuracy of the numerical approximation itself, 
which is sufficient for many applications. In the case of a problem with an incomplete 
set of end conditions, such as an initial value problem with condition e(a;(-l),a;(l)) = 
a;(—1) — x^ = D, K may be defined as 

I = 

= ■* 0, i = N, 

0, * 7^ 0, N. 

Lastly, for more complex BCs such as periodic conditions e{x(-l),x(l)) = a;(-l) - 
a;(l) = 0, or other complicated BCs such as nonlinear functions of x(—l) and a;(l), k 
may be defined as, Hi = 0, for i = 0,1,..., N, and the condition e(a;(-l),a;(l)) = 0 
may be enforced as a set of constraints by the NLP. However, this last case will result in 
strong enforcement of the BCs and some of the advantages of using the weak boundary 
formulation may be lost. 

The dynamical constraint becomes 

N 

'^DijX^^ + Ki-c^^ <S^, i = 0,l,..., N. 

oo 

The path constraints are approximated by 

h{x^\u^^) <S^ -1, i = 0,l,...,N. 
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Lastly, the cost functional J[x{-),u{-)] is approximated by the LGL quadrature rule, 

N 

i=0 

where x^ = x ^^,..., x^^~\ and ..., u^^~\. To allow for a practi¬ 

cal search area for the optimal solution the following constraints are added 

{x^^ G X, G [/, i = 0,1,..., N}, 

where X and U are the search regions that contain the optimal solution of the discretized 
nonlinear optimization. 

5.1.2. Computation Strategy 

In order to guarantee feasibility of the discretization. Theorems 5.1 and 5.2 show 
that a relaxation of the dynamical equality constraint to inequality is required, for GOCM- 
W and GOCM-W, respectively. 

5.I.2.I. Computation Strategy for GOCM-W 

The computational strategy of the GOCM-W is to find the feasible solution x^{t) G 
X and u^{t) G U for the following cases: 

Case 1. u{-) is piecewise and a;(-) G and piecewise C^, 

Case 2. u{-), i;(-) G and m>2, 

that minimizes 

J{x^{-),u^{-)) = j F{x^{t),u^{t))dt +E{x^{-l),x^{l)), (5.2) 
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subject to the Galerkin constraints 




(t) {x^{t) - f{x^{t),u^{t))) dt + Ki 


< MN- 




(5.3) 


\h^{x^{t),u^{t))\\^, < MN-^, 


where a = ^ and (m — 1), for Case 1 and 2, respectively; M is a constant independent of 
N-, 




{ h, h > 0, 
0, h <0; 


and 



a;^(—1) — x^, 
xf — 


i = 0, 

i = N, 


0 , 


where e(a;(—l), a;(l)) = x{—l)—x^,x{l)—X'l' ^=[0,0]^. 


(5.4) 


5.I.2.2. Computation Strategy for GOCM-W 

The computational strategy of the GOCM-W is to find the feasible solution x^ eX 
and E U for the following cases: 

Case 1. u{-) is piecewise and x{-) E and piecewise C^, 

Case 2. u{-), x{-) E m > 2 and i;(™“^)(t)is of bounded variation in t G 

[-1,1], 

that minimizes 


N 


u") 




x^\u^^)wi + E{x 


NO ;j.NN 
Jb 


(5.5) 


i=0 


subject to the Galerkin constraints 

II N 


'^^DijX^^ + Ki — 
j=0 


< MN-^, i 

OO 


h{x^\u^^) < MN-^ ■ 1 , 


= o,i,...,iv, 

i = 0,1,..., iV, 


(5.6) 
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where a = ^ and (m — 1), for Case 1 and 2, respectively; M is a constant independent of 
N-, 

x-^ — , i = N, (5.7) 

0, i7^0,iV, 

\ 

where e(a;(—1), a;(l)) = [a;(—1) — x^, a;(l) — a;7] ^ = [0,0]"^. 

5.1.3. Feasibility of Solutions 

In order to guarantee feasibility of the discretization, Theorems 5.1 and 5.2 intro¬ 
duced in [83] show that a relaxation of the dynamical equality constraint to inequality is 
required, for GOCM-W and GOCM-W, respectively. 

5.I.3.I. Feasibility of GOCM-W 

Theorem 5.1 (Feasibility of GOCM-W). Given any feasible solution 1 1 -)- (x, u), for Prob¬ 
lem B, consider the following two cases: 

Case 1. uf) is piecewise and x(-) G and piecewise C^, 

Case 2. «(•), xf) G 77™“^ and m > 2. 

Then, there exists a positive integer Nq such that, for any N > Nq, GOCM-W has a 
polynomial feasible solution, {x^ (t), (t)) such that 

\\x{t) -x^{t)\\^, < MN-^ and \\u{t) - < MN-^, 

where a = ^ and {m — 1), for Case 1 and 2, respectively; and M is a positive constant 
independent of N. 

Proof Let p{t) be the (iV — l)-th order truncated Legendre polynomial approximation of 
x{t). By Lemmas 4.1 and 4.2 there is a constant Cq independent of N, for any N > Ni, 
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such that 


\x{t) -p{t )\\^2 < CoN 


where a = | and (m — 1), for Case 1 and 2, respectively. Define 


x^(t)= / p{s)ds + x{—l). 


Thenp(t) = x^{t) and 


x{t) - x^{t)\\^^ < 2coN ", 


since, from Holder’s inequality (Lemma 4.3), we have 


x(t)—x^(t) =\ {x{s) — p{s)) ds < / |i;(s) — p(s)|(is 


<a/ 2 ( / |i;(s)-p(s)p(is ) = V2\\x{t) - p{t )\\^2 < V2coN 


Let u^{t) be the iV-th order Legendre polynomial so that 


u{t) - u^{t) 2 < CiN~ 


Recall that the LGL quadrature weights, {wk}k=o^ the property [46], 


Wk<C2N ^(1-( 4 )^) 2 , A; = 1,2,... ,iV - 1, 


for constant C2 > 0 independent of k and iV; for A; = 0, iV, 


Wo,Wn = 


N{N + l)' 
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From Theorem 4.1 we have, for each i = 1,..., iV — 1, 

J dt + Ki 

= J dt < CswfN-^^, 

for all N > N 2 , where C 3 is a constant independent of N. For i = 0, we have 

J 4>o W dt + Ko 

< j^(Poit) {x^{t) - f{x^{t),u^{t))) dt + |a;^(-l)-a;(-l)| 

< j 0^(f) < C3w|iV"“, 

since |a;'^(—1) — a;(—1) | = 0. For i = N,we have 

j {x^{t) - f{x^{t),u^{t))) dt + kn 

< j f{x^{i)^u^{i)))dt +\x^{l)-x{l)\ 

< cswiN-^ + V2coN-^, 

since 

|a;^(-l) -a;(-l)| < V2\\x{t) -pit )\\^2 < V2coN~°‘. 

Finally, for each i = 0,1,..., iV, 

J 4>k (t) ~ f (t)^(t))) dt + Ki < 

for all N > No, where M is a constant independent of N. 

The estimates for the path constraint follow from the proof of Theorem 4.1. 
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Thus a solution {t), {t)) to GOCM-W is feasible! 


□ 


5.I.3.2. Feasibility of GOCM-W 

Theorem 5.2 (Feasibility of GOCM-W). Given any feasible solution t h-)- (x, u ), for Prob¬ 
lem B, consider the following two cases: 

Case 1. uf) is piecewise and xf) G and piecewise C^, 

Case 2. m(-), x {-) G and m > 2. 

Then, there exists a positive integer TVq such that, for any N > Nq, GOCM-W has a feasible 
solution, {x^,u^) such that 


\\x{t) - x^{t)\\^, < MN-^, 

where a = ^ and {m — 1), for Case 1 and 2, respectively; and M is a positive constant 
independent of N. Additionally, u^(ti) = u{ti),fori = 0,1,..., iV. 

Proof Let p(t)be the {N — l)-th order truncated Legendre polynomial approximation of 
x{t) in the L^-norm. By Lemma 4.1 there is a constant do independent of N, for any 
N > Ni, such that 


\\x{t) -p(f)||icx> < doN 

where 13= {fn — for Case 2. For Case 1 we refer to [77-79] which show the truncated 
Legendre approximation for discontinuous functions with jump discontinuity (such as the 
step function defined in Lemma 4.2) displays Gibbs phenomenon. However the maximum 
amplitude of the overshoot has a finite limit; we conclude that for Case 1, = 0. Also, by 

Lemma 4.2 there is a constant di independent of N, for any N > N 2 , such that 

II^W -pWIIl2 < diN-^, 
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where a = ^ and (m — 1), for Case 1 and 2, respectively. Define 


if) = J p{s)ds + x{—l). 


Then p{t) = x^{t) and 


x{t) — x^(t) ^2 < 2(iiiV ", 


since, from Holder’s inequality (Lemma 4.3), 


x{t) — x^ (t) = 


(i;(s) — p{s)) ds 


'-1 


< / |i;(s) — p(s)|(is 




<a/ 2 ( / |i;(s)-p(s)p(is ) = V2\\x{t) - p{t)\\L 2 < V2diN~ 


Also, let u^{t) be an interpolating function of u{t). 


N 


where is any set of continuous functions with the property 

therefore 

= u{tj). 

Since x^{t) is a A^-th order polynomial, we have 


N 


X 


j=0 




6ij, and 
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where A is the (iV + 1) x (iV + 1) Legendre PS differentiation matrix (2.57) and 


Recall that the LGL quadrature weights, have the property [46], 

Wi < d2N~\l - i = 1,2,... ,N - 1, 



N 

'^DqjX^^ + Kq- 
j=0 
N 

< + |a;^(-l) - a;(-l)| 

3=0 

< d3WoN~^, 
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since 1) — a;(—1) | = 0. For i = N, 

N 

+ i^n — 

j=0 
N 

< + 1 ^( 1 ) - 

j=0 

= d^w^N-^ + V2d2N-^. 


Finally, for each i = 0,1,..., iV, 

N 

+ Ki- < MiV-“, 
i=o 

for all N > No, where M is a constant independent of N. 

The estimates for the path constraint follow from the proof of Theorem 4.2. 

Thus a solution {x^, u^) to GOCM-W is feasible! □ 

5.2. Galerkin Optimal Control with Continuous Element-based Galerkin 


We now consider a continuous element-based Galerkin approach. 

5.2.1. Method for Approximation 

In this approximation to Problem B, the weak integral form of x — f{x,u) = 0 in 
each element, takes the form [44] 


dx^^'>^(t) 


- f{x^^^^{t),u^^^^(t))]dt = 0, 


where Q = U^i defines the total domain. The state trajectory, x(t), is approximated 
inside each element, by interpolating iV-th order Lagrange polynomials, 
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at the nodes the relationship 


N 

j=0 

for e = 1,2,..., Ng, where are the LGL nodes, mapped back to the 

physical space inside each element, fie- Also, let u^{t) be an interpolating function of 

{u^^}U 

N 

j=0 


where {^jJj {t)}f=o are any set of continuous functions with the property V '>"(*.) = s.,. 


Therefore for e = 1, 2,..., A'e and j = 0,1,..., A^, and similarly, 

^(e)Nj ^ The relationship between the physical time domain, t G [to^f/] = 

, and the computational space, ^ ^ [—1,1], is given by [44] 


.( 1 ) ANe) 
I'O , I'N 




2 

At(^) 



1 and 


2 


dt, 


and conversely, 

t = (^ + 1) + Aq’ and dt = -^d^, 

where is the size of each element, fie, which can be nonuniform in length. 

The Lagrange polynomial defined on the LGL computational domain is given by 


N 

A (()=n 


i=o 

j¥=i 


(e- 0 ) 

te-0)’ 


f = 0 ,..., A. 
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The state trajectory, x, can now be approximated inside each element, by 


N 

j=0 

where ^re the Lagrange polynomials defined on the LGL grid. Likewise, u^{^) 

is given by 

N 

i=o 

where (6) = 

Remark 5.3. In this formulation = ^(e+i)^o ^{e)NN _ ^{e+i)m^ ^ _ 

1, 2,..., iVe — 1. This continuity condition is a consequence of the global formulation of 
the problem discussed in Remark 5.4. 


In the computational domain, the system becomes 

for e = 1, 2,..., iVe and i = 0,1,..., N, and in terms of the approximating polynomials 
becomes 


N 

j=0 ' 


kN 


'-1 


d^ ^ 2 




'-1 


In matrix-vector notation, our system can be expressed as 


N 

- 4 '' 

j=0 


= 0, i = 0,l,...,N, 


for e = 1, 2,..., iVe and f = 0,1,..., iV, where the local element (iV+1) x (iV+1) Galerkin 
differentiation matrix, is the same as that defined in Equation (4.3). If Q = iV LGL 
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quadrature nodes are used, the approximation to the (TV + 1) x 1 RHS vector simplifies to 



-{e)Ni 






0 , 1 , 


,iV, 


for e = 1,2,..., iVe and i = 0,1,..., TV, where the size of each element, can be 

nonuniform in length. 

Remark 5.4. So far the required objects have been identified to solve the system numer¬ 
ically with element-based Galerkin. However, since nodal basis functions are continuous 
across element boundaries and LGL nodes include both endpoints, a global solution to our 
problem can be found. To do this, a global assembly or direct stiffness summation can be 
done, where the direct stiffness summation operator is A^i- [44] 

The global equations to the problem become 

Np 

-c^^^ = 0, 1 = 1,..., Np. 

j=i 

The global Galerkin differentiation matrix, Djj and RHS vector, c^p^ are then defined by 

Djj = /\Dlf, and c^p^ = 

e=l e=l 


where Np = {NeN + 1) is the total number of grid points. Note that the direct stiffness 
summation operator does the mapping {i,e), (j, e) I, J [44]. See Section 2.3.2.1 for 
additional details. The dynamical constraint becomes 




J^Djjx^p^-c^ 


. 7=1 


<5, 1 = 1,2, 


,iV„ 
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The endpoint conditions and path constraints are approximated by 




oo 

<S-1, I = l,2,...,Np. 


< 6 


Lastly, the cost functional J[x{-),u{-)] is approximated by the LGL quadrature rule, 


J[a;(-), m(-)] 

A+(e) ^ 




e=l 


2=0 


where x^^ = ..., and u^p = [u^p^,u^p ‘^,..., u^p^p~\. To allow for a 

practical search area for the optimal solution the following constraints are included: x^p G 
X and u^p G U, where X and U are the search regions that contain the optimal solution 
of the discretized nonlinear optimization. 


5.2.2. Computation Strategy 

The computational strategy of the GOCM-CG is to find the feasible solution x^p G 
X and u^p G U that minimizes 

J^{x^p,u^p) 

e=l i=0 

+E{x^p^,x^p^p), 


subject to the Galerkin constraints 




y^Dijx 

i=o 


NpJ _ -NI 


\e{x p ,x 

I ^ ^ 11 OO 

h{x^p^ ,u^p^) 


<5^, 1 = 1,2,...,iVp, 

<5^-1, I = l,2,...,Np, 
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where 5^ is the feasibility tolerance, which is dependent upon N. 


5.3. GOCM with Discontinuous Element-based Galerkin 

We now consider a discontinuous element-based Galerkin approach introduced in 

[74], 


5.3.1. Method for Approximation 

In this approximation to Problem B, the weak integral form of x — f{x,u) = 0 in 
each element, yields [44] 

^ dt = 0 , 

where Vt = IJ^i defines the total domain. The state trajectory, x{t), is approximated 
inside each element, fie, by interpolating iV-th order Lagrange polynomials, 
at the nodes by the relationship 

N 

j=0 

for e = 1,2,, iVe, where are the LGL nodes, mapped back to the 

physical space inside each element, fie- Also, let u^{t) be an interpolating function of 

N 

j=0 

where are any set of continuous functions (not necessarily polynomials) with 

the property {ti) = <5^. Therefore for e = l,2,...,iVe and 

j = 0,1,..., iV, and similarly, The relationship between the physical 
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time domain, t G [to, t/] = tg^^ t^*"^ , and the computational space, ^ G [—1,1], is given 
by [44] 


— 1 and df = , , ^ dt, 


and conversely, 

t = (^ + 1) + and dt = -^d^, 

where At*^®) = t^^ — tg*^^ is the size of each element, which can be nonuniform in length. 
The Lagrange polynomial defined on the LGL computational domain is given by 

j¥=i 

The state trajectory, x, can now be approximated inside each element, fie, by 

N 

j=0 

where are the Lagrange polynomials defined on the LGL grid. Likewise, u^{^) 

is given by 

N 

i=o 

where = dij. In the computational domain, the system becomes 
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for e = 1, 2,..., iVe and i = 0,1,..., iV. Integration by parts on the first term yields the 
weak form relationship 



dcf). 


N 


di 


dS, + [ct>^ 


2 


i: 




0 . 


and in terms of our approximating polynomials, we have 


N 

J I d^ 


1 ^ , AA(e) /-I 

4>fdi +y; i4,r4'f];'4’ - ^ I = o. 


j=0 


'-1 


Remark 5.5. With the discontinuous element-based Galerkin approach, we let x, u and the 
basis functions be discontinuous across element edges. A numerical flux term x^*'^ acts as 
a jump condition between elements [44], Here, we consider the centered flux relationship, 
= I proposed by Delfour et al. [66], where e and q denote the element 

and its neighbor, respectively. 


Integrating by parts, yet again, results in the Galerkin strong form relationship 

^ /"l rlrh^ A+4) r^ 

j=o d^ 4-1 


for e = 1, 2,..., iVe and i = 0,1,..., iV. Since LGL nodes are used, the boundary term, 
may be simplified as 



{ I (a;(2)^o - xA)^^) , i = N, 

[o, i^N, 



1 (^^{Ne)m _ ^{Ne- 

0 , 


-l)ArAr^ 


i = 0, 
i 7 ^ 0 , 
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for elements fie = f^i and flAr^, respectively, and for each other element (fie 7^ fli, flArJ 
we have 


(e) 

Vi = <, 




X'- 




-l)iV7V^ 

)iV7V^ 


(e 


0 , 


i = 0, 
i = N, 
i^0,N. 


Remark 5.6. The problem’s endpoint conditions have not been introduced into the bound¬ 
ary term, and will instead be enforce in a strong sense through a set of endpoint 
constraints, as done in GOCM-S. If instead BCs are to be enforced weakly, a modification 
can be made to the boundary condition term, 


In matrix-vector notation, our system may be expressed as 


N 

- cf ^ = 0, e = 1, 2,..., iVe, * = 0,1,..., iV, 
i=o 


where the local element (iV -f 1) x (iV -f 1) Galerkin differentiation matrix, is the 
same as that defined in Equation (4.3). If Q = iV LGL quadrature nodes are used, the 
approximation to the (TV -f 1) x 1 RHS vector simplifies to 



-(e)Ni 






0,1,...,iV, 


where the size of each element, can be nonuniform in length. The dynamical con¬ 
straint becomes 


N 

j=0 






-{e)Ni 


<6, e = 1,2,... ,iVe, i = 0,1,... ,iV. 
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The endpoint conditions and path constraints are approximated by 




<S, 


< <^ ■ 1, e = 1, 2,...,iVe, i = 0,1,..., iV, 


Lastly, the cost functional J[x{-),u{-)] is approximated by LGL quadrature rule, 

e=l k=0 


where 




^iNe)m 

-(iVe)JViVj ^ 

U^= 





To allow for a practical search area for the optimal solution the following constraints are 
included: E X and E U, where X and U are the search regions that contain the 

optimal solution of the discretized nonlinear optimization. 


5.3.2. Computation Strategy 

The computational strategy of the GOCM-DG is to find the feasible solution x^ E 
X and E U that minimizes 




e=l 


i=0 
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subject to the Galerkin constraints 


N 




(e)^(e)Afj _ ^(e) _ -(e)Ni 


Vi - c' 


j=0 


<S^, e = l,2,...,Ne,i = 0,l,...,N, 


|g(5.(l)iVO^-(iV.)iV7V)||^<^iV^ 

■ 1, e = 1, 2,..., iVe, i = 0,1,..., iV, 


where 5^ is the feasibility tolerance, which is dependent upon N. 

5.4. Galerkin Optimal Control for Multi-scale Problems 

In this section, a multi-scale Galerkin optimal control approach is proposed to solve 
a specific optimal control problem, one in which the system dynamics are of different 
timescales (see Problem B). Such a problem may consists of a fast state, Xf{t), associated 
with the fast dynamics and a slow state, Xs{t), associated with the slow dynamics. GOCM- 
S (see Chapter 4) serves as the basis for the Galerkin multi-scale approach to Problem B. 
In particular, system (4.6) is fundamental to this alternative Galerkin optimal control for¬ 
mulation described below. The outline shown here follows the strategy given by Gong et 
al. in [71]. 

5.4.1. Method for Approximation 

The states and controls are approximated with globally interpolating Lagrange poly¬ 
nomials on different LGL timescales. The slow state, Xs{t), is approximated on sparse grid 
while the fast state, Xf{t), on dense grid where M < N. The slow and 

fast states are defined by the following approximating polynomials 

M 

Xsit) ^ xfit) = {t)xf^, 

j=0 
N 

Xf{t) ^ xf {t) = 

j=0 
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where the Lagrange polynomials {0^(t)}^Q and {0^(t)}^Q are defined on grids 
and respectively. Let 


xf^^XsiTj), j = 0, 

xf^^Xfitj), j = 0,1,...,N, 

and similarly, ^ u{tj), for j = 0,1,..., iV. 

Remark 5.7. For simplicity, the control variable, u(t), is approximated on the dense grid 
{ti}^Q, however this need not be the case. Modifications may be made to the process 
outlined in this section to cast the control onto a unique grid, such as sparse grid 
where M < N [71]. 


For GOCM-MS, a solution to the differential equations Xs = f{xs,Xf,u) and = 
g{xs,Xf,u) may be approximated by discretizing the slow dynamics over the dense grid 
with the following formulation 


M 

- cf 

j=o 


N 

E 

i=o 




-Ni 


0, f = 0, 

0 , f = 0,1 ,..., iV. 


The (TV + 1) X (M + 1), non-square, differentiation transformation matrix and (iV + 
1) X (iV + 1), square, differentiation matrix are defined by 


j-,NM _ 


j-)NN _ 




dt 


d(h^ 

dt = f = 0,1,..., iV, j = 0,1,. 




dt 

d(j)f 


M, 




{ti)wi = AfAwi, i,j = 0,1, 


TV, 
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where {wj}^Q are the LGL weights associated with LGL points is the 

standard {N + 1) x (N + 1) Legendre PS differentiation matrix (2.57) and is the 
{N + 1) X (M + 1) Legendre PS differentiation transformation matrix (2.39). The RHS 
vectors, cf and , are defined by 


cf = /(xf, xf\ = 0, * = 0,1,..., iV, 

= 0, i = 0,1,..., iV. 


The slow state approximation projected to the dense grid, may be calculated by the 
linear mapping with the relationship 


n 

S / J IJ s 

j=0 


for i = 0,1,..., iV, where is the {N + 1) x (M + 1) transformation matrix (2.37). 

Remark 5.8. Projecting the slow dynamics onto the dense grid provides a way of capturing 
the high frequency information of the fast state. If instead the intuitive approach is used, of 
discretizing the slow dynamics over the sparse grid, the high frequency information of the 
fast state is lost, resulting in a decrease in method accuracy [71]. 

The dynamical constraints therefore become 


M 



-Ni 

j=0 


N 



-Ni 

j=0 



i = 0,1,..., iV, 

i = 0,1,..., iV. 


The endpoint conditions and path constraints are approximated similarly by 


-MM -m -NN 


) <<^ 

' rv~i 


N 


h{xf\ xf, ■ 1, i = 0,1,..., iV. 
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Lastly, the cost functional J[x{-),u{-)] is approximated by the LGL quadrature rule, 


N 




.F{xf,xf,u^) = 


x^\ u^'')wi + E{x 

i=0 


-MO -MM -NO -NN\ 
s f s ) f 1 f )i 


where 






-MO -Ml 


^MM- 
■'S J 5 


x^ 


= X 


NO -N1 


, Xf 


?i.NNl 

’^/ J 


and 




NNl 


To allow for a search area for the optimal solution the following constraints are included: 
x^ e x^ G X J and G U, where Xg, Xj and U are the search regions that 
contain the optimal solution of the discretized nonlinear optimization. 


5.4.2. Computation Strategy 

The computational strategy of the GOCM-MS is to find the feasible solution x^ G 
Xg,xf E X f and G U that minimizes 


N 


Nt-M -N -N 


r{x 


Xf,U ) = 


i=0 


X 


Ni ^Ni r-,Ni 


, X 


f 


U 


)w. 


+ E{x 


MO 
s > 


X 


MM -NO ^NN 


, X 


f ’^/ 


subject to the Galerkin constraints 


M 



-Ni 

j=0 


N 



-Ni 

j=0 





<5^, 




f = 0,1,..., iV, 

f = 0,1,..., iV, 


h{ 


-Ni -Ni\ 
'^S ") f ’ ^ 




where 5^ is the feasibility tolerance, which is dependent upon N. 
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5.5. Modifications to Galerkin Optimal Control 


Two modifications to the general Galerkin optimal control formulation will be dis¬ 
cussed in this section: over-integration of the RHS vector and the use of quadrature points 
other than LGL, such as LG and LGR points. 

5.5.1. Over-Integration of the RHS Vector 

Thus far, all the Galerkin optimal control formulations discussed have approxi¬ 
mated the RHS vector with inexact integration, and N + 1 quadrature points. It may, 
however, be advantageous to approximate the RHS vector with increased accuracy. This 
may be accomplished by over-integration of the RHS vector while approximating the in¬ 
tegral using LGL quadrature. As outlined in the GOCM-S, the state trajectory, x{t), is 
approximated by 

N 

x{t) ^ x^{t) = 

j=0 

Let 


x^^^x{tj), j = 0, 

and similarly, is the approximation of u{tj). In the Galerkin optimal control formu¬ 
lation, a solution to the differential equation x — f{x, u) = 0 may be approximated at the 
LGL nodes with the following weak integral formulation [44] 

~ dt = 0 , 

and in terms of the approximating polynomials becomes 

j=o -^-1 
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In matrix-vector form the system becomes 

N 

-Ci = 0, i = 0,1,... ,N, 

j=0 

where the (iV -f 1) x (iV -f 1) differentiation matrix D is defined by Equation (4.2) and 
calculated exactly by Equation (4.3). The (iV + 1) x 1 RHS vector c is defined as 

d = j ^4)^ {t)f{x^{t),u^{t))dt, i = 0,l,...,N, 

and can be approximated with EGE quadrature by the relationship 

Q 

Ci ^ = '^4>f{tk)f{x^{tk),u^{tk))wk, i = 0,l,...,N, 

k=0 

where {wk}k=o the EGE quadrature weights (2.49) associated with EGE points, {tk}k=o- 

Remark 5.9. Recall that for LGL quadrature rule, Q = N + 1 integration points will 
integrate the RHS vector exactly when f{x{t),u{t)) is linear in x(t) and u(t). In the case 
of a nonlinear function, f, the accuracy of the numerical integration of the RHS vector may 
also be improved ifQ = N + l integration points are used when f consists of one or more 
nonlinear terms. 

Using over-integration, the RHS vector approximation becomes 

7V-I-1 

^Nt ^ '^(f)f{tfif[x^{tk),u^{tk))wk, i = 0,l,...,N, 

where {wk}k=o the EGE quadrature weights associated with EGE points, {ffcj^^.The 
dynamical constraint therefore becomes 

N 

i = 0,l,...,N. 
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The endpoint conditions and path constraints are approximated similarly by 

h{x^\u^^) < 5^ ■ 1, i 

where 1 denotes [1,..., 1]^. Lastly, the cost functional J[x{-),u{-)] is approximated by the 
LGL quadrature rule, 

N 

i=0 

where x^ = x^^ ■ ■ ■ x^^~\ and ■ ■ ■ u^^~\. To allow for a practical 

search area for the optimal solution the following constraints are included: x^ G X and 

e U. 

5.5.2. Galerkin Optimal Control with LG and LGR/F-LGR Quadrature Points 

The Galerkin weak formulation with weak boundary condition enforcement also 
allows for consideration of quadrature points that do not include both endpoints, e.g., LG 
and LGR (or F-LGR) nodes. This is advantageous since LG and LGR quadrature rules may 
lead to increased accuracies when performing the RHS vector integration. Recall that LG 
quadrature rule integration is exact for polynomial integrands of degree less than or equal 
to 2N + 1, where the LG nodes, defined by —1 < to < ■ ■ ■ < < 1 and 

are the roots of Equation (2.41). LGR (and F-LGR) quadrature rule is exact for polynomial 
integrands of degree less than or equal to 2N, where the LGR nodes, are defined 

by to = -1 < < ■ ■ ■ < tTV < 1, and are the roots of Equation (2.43) and the F-EGR 

nodes are the negative of the EGR points. 
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5.5.2.1. Method for Approximation for Galerkin Optimal Control with LGR/F-LGR 
Nodes 

In this section, a formulation is proposed to solve a specific optimal control prob¬ 
lem, one in which the initial or final conditions of the problem dynamics are provided (not 
both). Consider Problem B such that one of the two following cases exist: 

Case 1: e(a;(—1),a;(l)) = x{—l) — = 0, 

Case 2: e(a;(—1), a;(l)) = a;(l) — = 0, 

where x^ and x-^ are constants. As with the GOCM-W, we will consider the weak enforce¬ 
ment of the end conditions, however, now we will use the Galerkin weak form. In this 
approximation to Problem B, the state trajectory, x{t), is approximated with globally inter¬ 
polating 7V-th order Lagrange polynomials, defined on a grid of F-LGR or LGR 

nodes, Case 1 and 2, respectively, 

N 

x{t) ^ x^{t) = 

j=0 

Due to the property of the Lagrange polynomials, 4>^{U) = Sij, we have 

Also, let u^{t) be an interpolating function of 

N 

«"(i) = (ijfl"'. 

i-o 

where is any set of continuous functions with the property ipf {ti) = 6ij. Taking 

the weak integral form of x — f{x, u) = 0 yields [44] 

- fix-it)dt = 0, 
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for i = 0,1,..., iV. Integration by parts on the first term results in Galerkin weak form, 


dfj). 


'-1 


dt 


dt + [ct>f 


- (l)^f{x^,u^)dt = 0. 


'-1 


In terms of our approximating polynomials, we have 




for f = 0,1,..., iV. The expression may be simplified as 


'^DijX^^ + Ki- c^* = 0, 
j=0 

for each f = 0,1,..., iV, where 

l)a;°, i ^ N, 

l)a;° + x^^, i = N, 

for Case 1 and 

j4 )q — x^^, 

for Case 2. 

The {N + 1) X {N + 1) differentiation matrix D is defined by 
Dii = - *, J = 0,1,..., /V. 


i = 0, 
i 7^ 0, 
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If F-LGR/LGR quadrature rule is used with Q = N quadrature points, the differentiation 
matrix, D, ean be ealeulated exaetly by the relationship 


n — '^d(t)f{tk) 

~ ^ dt 

k=0 


{tk)Wk = 


d<Pi 


N 


dt 








i,j = 0, 


where {wj}f^Q are the F-LGR/LGR quadrature weights (2.50) and Aij = is the 

F-LGR/LGR PS differentiation matrix. 

The RHS veetor, c ean be approximated by the relationship 

Q 

a ^ = '^4>^{tk)f{x^{tk),u^{tk))wk = f{x^\u^")wi, i = 0,l,...,N. 

k=0 

Remark 5.10. Recall that for F-LGR/LGR quadrature rule, integration is exact for polyno¬ 
mial integrands of degree less than or equal to 2N. IfQ = N integration points are used, 
the RHS vector will integrate exactly when f {x(t), u(t)) is linear in x(t) and u(t). In the 
case of a nonlinear function f, the accuracy of integration may be improved by increasing 
the number of quadrature points Q (See Section 2.2.2). 

The dynamieal eonstraint beeomes 


N 

+ ki- 


j=0 


< S 


N 


/ = 0,1,..., iV. 


The path eonstraints are approximated by 


h{x^\u^^) <6^ -1, / = 0,l,...,iV. 


Lastly, the eost funetional .J[x{-),u{-)] is approximated by the F-LGR/LGR quadrature 
rule, 

N 

i=0 
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where ■ ■ ■ x^^~\ and ■ ■ ■ u^^~\. To allow for a practical 

search area for the optimal solution the following constraints are included: x^ G X and 

e U. 

Remark 5.11. Note that the use of the Galerkin optimal control formulation outlined here 
with LGR or F-LGR nodes does not automatically provide the optimal control solutions at 
one endpoint of the domain (this applies to t = 1 for LGR points and t = —1 for F-LGR 
points). However, if required by the application, the control solutions at t = 1 or t = —1 
may be found by interpolation of the control approximation, with a possible reduction 
in accuracy. 

5.5.2.2. Method for Approximation for Galerkin Optimal Control with LG Nodes 

In this section, a formulation is proposed to solve a specific optimal control prob¬ 
lem, one in which a complete set of boundary conditions are provided for the problem 
dynamics. Consider Problem B such that 

e(a;(— 1), a;(l)) = [a;(— 1) — x^,x{l) — x-^] = [0, 0]^ , 

where x^ and x-^ are constants. As with the GOCM-W, we will consider the weak enforce¬ 
ment of the boundary conditions, however, now we will use the Galerkin weak form. In 
this approximation to Problem B, the state trajectory, x{t), is approximated with globally 
interpolating 7V-th order Lagrange polynomials,{0^}^O’ defined on a grid of LG nodes, 

{t,}U 

N 

x{t) ^ x^{t) = {t)x^L 

j=0 

Due to the property of the Lagrange polynomials, 4>^(U) = Sij, we have 

x^^=x^{tj), j = 0,l,...,N. 
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Also, let u^{t) be an interpolating function of 

N 

U^it) = 

j=0 

where is any set of continuous functions (not necessarily polynomials) with the 

property = 5ij. Taking the weak integral form of i; — f{x, u) = 0 yields [44] 

- fix-it)dt = 0, 

for i = 0,1,..., A^. Integration by parts on the first term results in Galerkin weak form, 



In terms of our approximating polynomials (and introducing the true initial condition, 
a;^(—1) —)■ a;(—1) and a;'^(l) —)■ a;(l)) we have 

for i = 0,1,..., A^. By letting a;(—1) = and a;(l) = x^, the expression may be 
simplified as 

N 

'^^DijX^^ + ki — = 0 , 

j=0 

for each i = 0,1,..., A^, where 
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The (TV + 1) X (iV + 1) differentiation matrix, D, is defined by 


Dij = - 




dt 


i,j = 0,1,...,N. 


If LG quadrature rule is used with Q = N quadrature points, the differentiation matrix, D, 
can be calculated exactly by the relationship 


n — 

^ dt 

k=0 


{tk)Wk = 


dc^l 


N 


dt 








hJ = 0, l,...,iV, 


where the LG quadrature weights (2.48) and Aij = <j)f{ti) is the LG PS 

differentiation matrix. 

The RHS vector, c can be approximated by the relationship 

Q 

a ^ {tk)f{x^{tk),u^{tk))wk = f{x^\u^^)wi, i = 0,l,...,N. 

k=0 

Remark 5.12. Recall that for LG quadrature rule, integration is exact for polynomial in¬ 
tegrands of degree less than or equal to 2N -\-1. If Q = N integration points are used, the 
RHS vector will integrate exactly when f{x{t),u{t)) is linear in x(t) and u(t). In the case 
of a nonlinear function f, the accuracy of integration may be improved by increasing the 
number of quadrature points Q (See Section 2.2.2). 


The dynamical constraint becomes 


TV 


j=0 


< S 


N 


i = 0,1,..., N. 


The path constraints are approximated by 


h{x^\u^^) <6^ -1, i = 0,l,...,N. 
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Lastly, the cost functional J[x{-),u{-)] is approximated by the LG quadrature rule, 


N 

i=0 

where x^ = x^^ ■ ■ ■ x^^~\ and ■ ■ ■ u^^~\. To allow for a practical 

search area for the optimal solution the following constraints are included: x^ G X and 
e U. 

Remark 5.13. Note that the use of the Galerkin optimal control formulation outlined here 
with LG nodes does not automatically provide the optimal control solutions at the endpoints 
of the domain, t = ±1. However, if required by the application, the control solutions at 
t = —I and/or t = 1 may be found by interpolation of the control approximation, with 
a possible reduction in accuracy. 
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CHAPTER 6: 

GALERKIN OPTIMAL CONTROL WITH LEGENDRE 
POLYNOMIAL TEST EUNCTIONS 


In this chapter, a Galerkin optimal control formulation is proposed where Legen¬ 
dre polynomials replace Lagrange polynomials as test functions in the weak integral ap¬ 
proximation of the problem dynamics. The purpose of this modification is highlighted in 
consistency Theorem 6.2, which is valid for problems with discontinuous controls (unlike 
consistency Theorems 4.3 and 4.4 for Problems GOCM-S and GOCM-S, respectively). 
Theorem 6.2 proves that the nonlinear programming Problems GOCM-L is a consistent 
approximation to the continuous optimal control Problem B, even those with piecewise 
continuous controls. 

6.1. Methods for Approximation 

In the Galerkin optimal control approximation to Problem B, with Legendre poly¬ 
nomial test functions, the state trajectory, x{t), is approximated with globally interpolating 
A^-th order Lagrange polynomials, defined on a grid of LGL nodes, 

N 

x{t) ^ x^{t) = {t)x^K 

3=0 

Due to the property of the Lagrange polynomials, (j)^{U) = Sij, we have 

Also, let u^{t) be an interpolating function of 

N 

u^it) = 

3=0 
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where is any set of continuous functions (not necessarily polynomials) with the 

property = 5ij. In this formulation, a solution to the differential equation x — 

f{x,u) = 0 may be approximated at the LGL nodes with the following weak integral 
formulation 

for i = 0,1,. .., iV, where the test functions, Li = are the normalized Legendre 

polynomials of order i. The L^-norm of Lj is given by [72] as 

2i + l‘ 

In terms of the approximating polynomials, the system becomes 

/ Li—^dtx^^— / Lif{x^,u^)dt = 0. 

In matrix-vector form, the system becomes 

N 

= z = 0,l,...,N, 

j=0 

where the (iV + 1) x (iV + 1) differentiation matrix is defined by 

j. ~ d(j)f{t) 

and the (iV + 1) x 1 right-hand-side (RHS) vector c is defined as 

cf = [ Li{t)f{x^{t),u^{t))dt, i = 0,l,...,N. 
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If LGL quadrature rule is used, with Q = N quadrature points, the differentiation matrix, 
D^, can be calculated exactly by the relationship 


_ d©;.' 


If Q = N quadrature points are used, the RHS vector, c^, can be approximated by the 
relationship 


^ -Ni _ 


' =y^Li{tk)fix^^, u^^)wk, i = 0,1,..., iV. 


The system may be simplified as 


-c^^ = 0, z = 0,l,...,N. 


The dynamical constraint therefore becomes 


Si", » = o,i,... 


The endpoint conditions and path constraints are approximated similarly by 




h{x^\u^^) <6^ - 1 , f 


Lastly, the cost functional J[x{-),u{-)] is approximated by the LGL quadrature rule. 
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where ■ ■ ■ x^^~\ and ■ ■ ■ u^^~\. To allow for a practical 

search area for the optimal solution the following constraints are included: x^ G X and 
G U, where X and U are the search regions that contain the optimal solution of the 
discretized nonlinear optimization. 

6.2. Computation Strategy 

The computation strategy for Galerkin optimal control with Legendre polynomial 
test functions is presented in two forms. First, the strategy for the continuous problem, 
in terms of the approximating polynomials is outlined, denoted as GOCM-L. Next, the 
discrete problem, discretized on a LGL grid is presented, denoted as GOCM-L. 

6.2.1. Computation Strategy for GOCM-Z 

The computational strategy of the GOCM-L is to find the feasible solution x^{t) G 
X and u^{t) G U for the following cases: 

Case 1. u{-) is piecewise and x{-) G and piecewise C^, 

Case 2. u{-), x(-) G and m>2, 

that minimizes 

J{x^{-),u{-)) = j F {x^{t),u{t)) dt + E {x^{-l),x^{l)) , 

subject to the Galerkin constraints 

< 




Li{t) {x^{t) - f {x^{t),u{t))) dt 

\\e (a;^(-l),a;^(l)) 
||/i+ {x^{t),u^{t))\ 


L2 
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where a = ^ and (m — 1), for Case 1 and 2, respectively; M is a constant independent of 
N and 


{ h, h > 0, 

0 , h<0. 


6.2.2. Computation Strategy for GOCM-L 

The computational strategy of the GOCM-L is to find the feasible solution x^(t) E 
X and u^(t) E U that minimizes 


N 


u^) = {x^\ w, + E {x^^, x^^) , 


i=0 


subject to the Galerkin constraints 


N 




L^Nj _ 


j=0 


<S^, i = 0,l,...,N, 


|e(x-x--)|L<^-, 

h{x^\u^^) <S^ -1, i = 0,l,...,N. 


6.3. Feasibility of Solutions 

Theorem 6.1 (Feasibility of GOCM-Z). Given any feasible solution t h- )■ {x, u), for Prob¬ 
lem B, consider the following two cases: 

Case 1. uf) is piecewise and xf) E and piecewise C^, 

Case 2. «(•), xf) E and m > 2. 

Then, there exists a positive integer Nq such that, for any N > Nq, GOCM-L has a poly¬ 
nomial feasible solution, {x^(t), (t)) such that 


x{t)-x^{t)\\^, < MN-'^ and \\u{t) - {t)\\< MN-^, 
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where a = ^ and {m — 1), for Case 1 and 2, respectively; and M is a positive constant 
independent of N. 

Proof Let p{t) be the (iV — l)-th order truncated Legendre polynomial approximation of 
x{t). By Lemmas 4.1 and 4.2 there is a constant cq independent of N, for any N > Nq, 
such that 


-pWIIl2 < CqN 


where a = ^ and (m — 1), for Case 1 and 2, respectively. Define 


(t) = J p{s)ds + x{—l). 


Then p{t) = x^{t) and 


x{t) - x^{t )\\^2 < 2coiV ", 


since, from Holder’s inequality (Lemma 4.3), we have 


x{t) — x^ (f) = 


{x{s) — p{s)) ds 




< / |i;(s) — p(s)|(is 


'-1 


<a/2 ( / |i;(s)-p(s)|^(is ) = A/2||i;(t) - p(t)||^2 < v^coiV 


( 6 . 1 ) 


Let u^{t) be the iV-th order Legendre polynomial so that 


u{t) - M"(t) ^2 < CiN ". 
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From our Galerkin approximation, Holder’s inequality (Lemma 4.3), and the property, 

L2 = i = 0,1,..., iV, 


< 


r 

-1 

Ht) 


Li{t) {x^{t) - f {x^{t),u^{t))) dt 


2i+l 


L2 


= 11^2 
< ||i;(t) -i;^(t)||^2 + 1|/ {x{t),u{t)) - f{x^{t),u^{t))\\^^ 

= coiV“5 + /i||a;(t) - a;^(t )||^2 + l2\\u{t) - 

< CoAr-“ + 2/iCoiV-“ + /2CiiV-“, 

where li and I 2 are the Lipschitz constants of / with respect to x and u, respectively, which 
are independent of N. It follows that 




Li{t) {x^{t) - f {x^{t),u^{t))) dt 


< MN-^, 


and holds for each i = 0, 1, ..., iV, and all N > Nq, where M is a constant independent of 
N. 

For the endpoint condition we have 


a;(l) — = 


(i;(s) — p{s)) ds 


1-1 


< / |i;(s) — p(s)|(is 




<^2 ( ^ |i;(s)-p(s)p(is] = V2\\x(t) - p(t )\\^2 < V2coN , 


so we have, by Lipschitz condition. 
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For the path constraint let P = [t\h (t), (t)) > O} , = [—1,1] \T>, since 

h{x{t),u{t)) < 0. Then 


< 


/i+ {x^{t),u^{t))\\^2 = d?j 

< {h {x^{t)^ — h{x{t),u{t))Y 


{h{x^ (t)) — h{x{t),u{t))ydt + / {h{x^ (t)) — h{x{t),u{t))ydt 


'V 


IV 


= (^h (^x^{t),u^{ty — h{x{t),u{t))Y dt^ 

= \\h{x^{t),u^{t)) - h{x{t),u{t))\\^, 

< hWxit) -x^{t)\\^^ + k\\u{t) -M^(t)||^2 < MN~'^, 


1 

2 


where /s and /4 are the Lipschitz constants of h with respect to x and u, respectively, which 
are independent of N. Hence 


||/i+ <MiV-“ 

Thus a solution [x^{t), {t)) to GOCM-L is feasible! □ 

6.4. Consistency of Solutions 

Theorem 6.2 (Consistency of GOCM-L). Suppose [x^ (t), (t)) is a solution of GO CM- 

L and there exists {x(t),u(t)) such that: 

Case 1. u{-) is piecewise and xf) G and piecewise C^, 

Case 2. u{-), x{-) G and m > 2. 

Also, suppose 


lim \u{t) — (t) = 0) 

N^oo ' 


( 6 . 2 ) 
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and x^(t) —)■ x(t) uniformly, thus 


lim x{f) — x^{t) 2 = 0- 

N^OO ^ 


(6.3) 


Then {x{t)^u{t)) satisfies 

Mt) - f = 0, 

< e(a;(-l),a;(l)) = 0, 
h {x{t),u{t)) < 0, 

and is an optimal solution to Problem B. 


Proof. Due to the completeness of the Legendre polynomials [80], to prove \\x(t) — f {x, u) ||^2 
0 it is sufficient to prove 



Li{t) {x{t) - f {x{t),u{t))) dt 


0 , 


for each f = 0,1,..., cx) (see definition of completeness, Definition 4.2). Consider 


Li{t) {x{t) - f {x{t),u{t))) dt 


'-1 


< 


i: 


+ 




Lft) [x{t) — x^{t)) dt 


Lft) {x^{t) - f {x^{t),u^{t))) dt 

Li{t) (/ {x^{t),u^{t)) - f {x, u)) dt 
< MN-'^ + ||i;(t) -x^{t )\\^2 + \\f {x{t),u{t)) - f {x^{t),u{t ))\\^2 

= MN~^ + ||a;(t) - a;^(t)||^2 + (il|a^(^) - a^^(^)|L2 + h\\u{f) - M^(t)||^2, 




where a = \ and (m — 1), for Case 1 and 2, respectively (from Theorem 6.1); M is a 
positive constant and li and h are the Lipschitz constants of / with respect to x and u. 
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respectively, all independent of N. It follows that as iV —)■ cx) we have 


- f {x{t),u{t ))\\^2 = 0. 

For the endpoint condition, since x^{t) —)■ x{t) uniformly, we have —)■ a;(l) 

and a;^(—1) —)■ a;(—1). Since, from the formulation of the computational strategy we have 
|e 1), a;^(l)) | < we conclude that e (a;(—1), a;(l)) = 0 as iV —)■ cx). 

For the path constraint, since h{x{t),u{t)) is piecewise if h{x{t*),u{t*)) >0, 
3 an interval (a, h) in which h {x{t),u{t)) > 0. Then 

{h{x{t),u{t))f d?j >0. 

However, 

\\h{x{t),u{t))\\^ 2 ^^^f,j < \\h{x{t),u{t)) - h+ + \\h+ {x^(t), (t)) 

< \\h {x{t),u{t)) - h {x^(t),u^(t)) + MiV“ 

< hWxit) - a;^(t))||^2 + l4\\u{t) - M^(t))||^2 + MN~'^, 

where and I 4 are the Lipschitz constants of h with respect to x and u, respectively, 
which are independent of N. Hence, this is a contradiction, therefore h {x{t),u{t)) < 0 as 
N ^ 00 . 

Suppose that {x{t),u{t)) is not optimal. Then 3 {x*(t), u*(t)) so that 
J {x*{-),u*{-)) < J {x{-),u{-)). 

Also, 3 {x*(t), u*(t)) such that 

||a;*^(t) - a;(t )||^2 < and - u{t)\\^, < MN-^, 


L‘2{a,b) 
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where (t), u*^(t)) is a feasible trajectory of GOCM-L. Therefore 

(6.4) 

However, 

\J - J (a;(.),M(-))| 

< j \F {x^{t),u^{t)) - F {x{t),u{t))\dt + |E (a;^(-l),a;^(l)) - a;(-l),a;(l))| 

< V2\\F {x^{t),u^{t)) - F {x{t),u{t)) dt\\^^ + |E (a;^(-l),a;^(l)) - a;(-l), a;(l))|. 

Due to the Lipschitz condition and assumptions (6.2) and (6.3) we have 

- J (a^(-),«(■))! = 0. 

Similarly, 

- J (a^*(-),w*(-))| = 0. 

Therefore, from (6.4) we have 

J {x*{-),u*{-)) > J {x{-),u{-)). 

This is a contradiction, since we assumed 


J {x*{-),u*{-)) < J (a;(.), «(■)). 


We conclude that {x{t),u{t)) achieves an optimal cost and therefore is an optimal solution 
to Problem B! □ 
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Remark 6.1. Theorems 6.1 and 6.2 show that solutions to the GOCM-L exist and provide 
confidence that they will converge to the optimal solutions of Problem B. More importantly, 
Theorem 6.2 provides a foundation to show that solutions to the GOCM-L will converge to 
optimal solutions with discontinuities in the control, such as solutions to bang-bang control 
problems. However, as with the GOCM-S, questions still remain about the conditions 
under which the underlying assumptions exist (in the case of the GOCM-L, assumptions 
(6.2) and (6.3)), but the required analysis is above the scope of this dissertation. 
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CHAPTER 7: 

EXAMPLE PROBLEMS 


A number of examples are provided in this chapter to highlight the versatility and 
accuracy of the Galerkin optimal control formulations discussed in Chapters 4-6. Exam¬ 
ples 7.1, 7.2 and 7.4 demonstrate the improved accuracy provided by the Galerkin optimal 
control formulation with weak enforcement of boundary conditions over strong enforce¬ 
ment for problems with fixed boundary conditions as well as incomplete sets of end con¬ 
ditions. In particular, the examples show that the GOCM-W has an advantage over the 
GOCM-S for low order approximations of control solutions. Examples 7.1 and 7.4 also 
show the potential advantages of the Galerkin formulations with E-EGR and EG points, re¬ 
spectively. Additionally, Example 7.3 demonstrates the effectiveness of the element-based 
Galerkin optimal control formulations (such as the GOCM-DG) when employed to approx¬ 
imated optimal control problems with discontinuous controls. Easily, Examples 7.5 and 7.6 
demonstrate the computational efficiency in which the multi-scale Galerkin optimal control 
formulation may solve multi-scale problems, those with states and controls that evolve on 
different timescales. In contrast to the difficulties with the multi-scale Eegendre PS method 
(see Section 3.4) highlighted in Chapter 3, the GOCM-MS is shown to successfully reduce 
the size of multi-scale problems. 

7.1. Example 7.1: Nonlinear Two-Dimensional Problem with Eixed 
Initial Conditions 

Consider the nonlinear two-dimensional problem wiih. fixed initial conditions given 
by Gong et al. [3] of minimizing the cost function 

J = Axi(2) + X‘2{2) + A f v?dt, (7.1) 

Jo 
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subject to the dynamics 


Xiit) = x\{t) and X2{t) = u{t), 


and with initial conditions 


a;i( 0 ) = 0 and 0:2(0) = 1 . 


The analytic solution to this problem is given by 


Xiit) 

X2{t) 

u{t) 


2 64 

5 ~ 5(2 + t)5’ 
4 

(2TtF’ 

-8 

(2TtF’ 


obtained via Pontryagin’s maximum principle. This problem was solved using the GOCM- 
W (Section 5.1) with optimality and feasibility tolerances of 10“® and 10“®, respectively. 
A two-point initial guess was provided. Figure 16 shows a comparison of the exact solution 
with the GOCM-W approximation and N = 20. 
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Figure 16: Exact solution and GOCM-W approximation with = 20 for Example 7.1. 
Figures 17, 18 and 19 show the maximum error of the states and control, 


error^illoo = 



i = 0,1,. 


error^ailoo = 

\x2{ti) 


i = 0,1,. 

..,N, 

||error„||^ = | 

\u{ti) - 


i = 0,1,.. 

.,N, 


respectively, for the LPM (Legendre PS method. Section 3.2), GOCM-S (Chapter 4), 
GOCM-W and GOCM-FLGR (Section 5.5.2.1) vs. polynomial order, N. Note that the 
GOCM-S approximations show nearly the exact same exponential convergence rates as the 
Legendre PS method numerical solutions throughout the displayed range of N values. 
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Figure 17: State x\ approximation error vs. polynomial order, N, for Example 7.1. 


However, the GOCM-W and GOCM-FLGR numerical solutions of X 2 and u show a marked 
improvement over that of the GOCM-S. Observations from this example show that the 
GOCM-W and GOCM-FLGR formulations have the potential to provide increased accura¬ 
cies of the control solution with lower order approximations. 


162 









Figure 18: State X 2 approximation error vs. polynomial order, N, for Example 7.1. 



Figure 19: Control u approximation error vs. polynomial order, N, for Example 7.1. 
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7.2. Example 7.2: Nonlinear Two-Dimensional Problem with Fixed 
Boundary Conditions 


Consider the nonlinear two-dimensional problem with fixed boundary conditions 
given by Kang [6] of minimizing the cost function 


J = / (1 — Xi -f X 1 X 2 + Xiu^dt, (7.2) 

Jo 


subject to the dynamics 

Xi(t) = —x\x 2 and X 2 (t) = —1 -\ - 

and with boundary conditions 

a;i(0) = 1, 0 : 2 ( 0 ) = 0, a:i(7r) = 


+ X2 + sin t + u, 


and 0:2 (tt) = 2. 


The analytic solution to this problem is given by 


Xiit) 

X2{t) 

u{t) 


1 

1 — sin t t ’ 

1 — cos t, 

— (t + l)+ sin t + cos t, 


obtained via Pontryagin’s maximum principle. This problem was solved using the GOCM- 
W (Section 5.1) with optimality and feasibility tolerances of 10“® and 10“®, respectively. 
A two-point initial guess was provided. Figure 20 shows a comparison of the exact solution 
with the GOCM-W approximation and N = 20. 
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Figure 20: Exact solution and GOCM-W approximation with = 20 for Example 7.2. 
Figures 21, 22 and 23 show the maximum error of the states and control, 


error^illoo = 



i = 0,1,. 


error^ailoo = 

\x2{ti) 


i = 0,1,. 

..,N, 

||error„||^ = | 

\u{ti) - 


i = 0,1,.. 

.,N, 


respectively, for the LPM (Legendre PS method. Section 3.2), GOCM-S (Chapter 4), 
GOCM-W and GOCM-OI (Section 5.5.1) vs. polynomial order, N. Note that the GOCM- 
W approximation represents a formulation for which the initial conditions are enforced 
weakly and the final conditions are imposed in a strong sense through an endpoint con¬ 
straint. Although enforcing all boundary conditions weakly provides accurate solutions, 
a partial weak enforcement of end conditions produces higher accuracies. Also note that 
the GOCM-OI approximations represent the over integration of this GOCM-W formula¬ 
tion. The GOCM-S approximations show nearly the exact same exponential convergence 
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rates as the Legendre PS method numerical solutions throughout the displayed range of N 
values. 



Figure 21: State xi approximation error vs. polynomial order, N, for Example 7.2. 

The GOCM-W numerical solutions of X 2 and u show a marked improvement over that 
of the GOCM-S for N < 22; the control error for the GOCM-W is up to two orders of 
magnitude lower than that of the GOCM-S for the N values in this range. Finally, the 
GOCM-OI formulation has a smoothing effect on the accuracies of the GOCM-W approxi¬ 
mations. While the GOCM-OI approximations of X 2 appear to be less superior to that of the 
GOCM-W, the GOCM-OI approximations of u gain slightly in accuracy for the lower order 
approximations. Observations from this example show that the GOCM-W and GOCM-OI 
formulations have the potential to provide increased accuracies of the control solution with 
lower order approximations. 
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Figure 22: State X 2 approximation error vs. polynomial order, N, for Example 7.2. 



Figure 23: Control u approximation error vs. polynomial order, N, for Example 7.2. 
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7.3. Example 7.3: Two-Dimensional Bang-Bang Control Problem with 
Fixed Boundary Conditions 

Consider the 2-dimensional bang-bang control problem with fixed initial conditions 
given by Pinch [21] of controlling the system from the initial point (a, &) at t = 0 to the 
origin (0, 0) in the shortest time. Here, we minimizing the cost function 

C/ 

J = / dt = tf, 

Jo 

subject to the dynamics 


x(t) = y(t) and y(t) = u(t), 

and with boundary conditions 

a;(0) = a, y{0) = b, x{tf) = 0, y{tf) = 0, 

and control constraint \u{t)\ < k. For the case where a = 1, & = 3 and k = 1, the analytic 
solution to this problem is given by 


x{t) = 


y{t) = 


u{t) = 


— y + 3f-|-l, t ^ 


—t -f 3, t < 

iy ^ 5 

- 1 , t<t^, 

1, t > 


obtained via Pontryagin’s maximum principle, where = 3 -f y y and tf = ?> + 2^ 
This problem was solved using the GOCM-S (Chapter 4), GOCM-W (Section 5.1), GOCM- 
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CG (Section 5.2) and GOCM-DG (Section 5.3). A two-point initial guess was provided for 
each approximation. 

Figures 24 and 25 show comparisons of the exact solution with the GOCM-S and 
GOCM-W numerical solutions, respectively, with approximation order, = 20. 



Figure 24: Exact solution and GOCM-S approximation with = 20 for Example 7.3. 

Erom observations, we can see a slight improvement from the GOCM-S to the GOCM- 
W numerical solutions, particularly in the approximation of the control, u, in the vicinity 
of the discontinuity location, However, the GOCM-S and GOCM-W have difficulty 
in approximating the discontinuous control solution. Even with an increased approxima¬ 
tion order, there remains a maximum error of 0(10“^) for the GOCM-S and GOCM-W 
approximations of the control. 
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Figure 25: Exact solution and GOCM-W approximation with = 20 for Example 7.3. 

Eigures 26 and 27 show comparisons of the exact solution with the GOCM-CG and 
GOCM-DG approximations, respectively, with number of nonuniform elements, N^. = 2, 
polynomial order inside each element, N = 10, and the boundary of the two elements 
located at It should be noted that in the element based formulations, is defined as 
a decision variable in the NEP. An initial guess for is then provided to the NEP with 
bounds prescribed. The total number of points for the GOCM-CG approximation is, Np = 
{Nf,N + 1) = 21 while for the GOCM-DG approximation, Np = {N + 1) Ae = 22. 
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Figure 26: Exact solution and GOCM-CG approximation with Np = 21, = 2 and the 

boundary of the elements located at for Example 7.3. 


Note that both the GOCM-CG and GOCM-DG approximations of the states, x and y, 
achieve maximum errors of 0(10“®) from the exact states. Additionally, the GOCM-DG 
achieves an impressive 0(10“®) maximum error from the exact control, u, as compared 
with an accuracy of 0(10“^) for the GOCM-CG. 
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Figure 27: Exact solution and GOCM-DG approximation with Np = 22, = 2 and the 

boundary of the elements located at for Example 7.3. 

7.4. Example 7.4: Two-Dimensional Problem with Fixed Boundary 
Conditions 

Consider again the linear two-dimensional problem wi\h fixed boundary conditions 
given in Example 3.1 of minimizing the cost function 



(7.3) 


subject to the dynamics 


xiit) = X 2 and X 2 {t) = C sm{kt) + u, 


and with boundary conditions 


a;i(0) = 0, 0 : 2 ( 0 ) = 0, xi{tf) = 1 and X 2 {tf) = 0. 
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The analytic solution to this problem is given by Equations (3.9)-(3.13), where tf = 10, 
C = 0.1 and k = 8. This problem was solved using the GOCM-W (Section 5.1) with 
optimality and feasibility tolerances of 10“® and 10“®, respectively. A two-point initial 
guess was provided for each approximation. 



Figure 28: Exact solution and GOCM-W approximation with = 50 for Example 7.4. 


error^illoo = 

\xiiti) 


i = 0,1,. 

..,N, 

error^ailoo = 

\x2(ti) 


i = 0,1,. 

..,N, 

||error„||^ = 

\u{ti) - 


i = 0,1,.. 

.,N, 


respectively, for the EPM (Eegendre PS method. Section 3.2), GOCM-S (Chapter 4), 
GOCM-W and GOCM-EG (Section 5.5.2.2) vs. polynomial order, N. Note that, as with 
Example 12, the GOCM-W approximation represents a formulation for which the ini¬ 
tial conditions are enforced weakly and the final conditions are imposed in a strong sense 
through an endpoint constraint. Although enforcing all boundary conditions weakly pro- 
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vides accurate solutions, a partial weak enforcement of end eonditions produces higher 
aeeuraeies. Note that the GOCM-S approximations show nearly the exaet same exponen¬ 
tial eonvergenee rates as the LPM numerieal solutions throughout the displayed range of 
N values. 



Figure 29: State xi approximation error vs. polynomial order, N, for Example 7.4. 
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Figure 30: State X 2 approximation error vs. polynomial order, N, for Example 7.4. 



Figure 31: Control u approximation error vs. polynomial order, N, for Example 7.4. 
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However, the GOCM-W and GOCM-LG numerical solutions for u show a marked im¬ 
provement over that of the GOCM-S, particularly in the lower order approximations; the 
control error for the GOCM-LG is one to five orders of magnitude lower than that of the 
GOCM-S in this range. Observations from this example show that the GOCM-W and 
GOCM-LG formulations have the potential to provide increased accuracies of the control 
solution with lower order approximations. 

7.5. Example 7.5: Two-Dimensional Multi-scale Problem with Fixed 
Boundary Conditions 

Consider again the linear two-dimensional problem with, fixed boundary conditions 
given in Example 3.1 and 7.4. Here we take advantage of the multi-scale nature of the 
problem. The slow state, xi, fast state, X 2 and control u are discretized on LGL grids, 

respectively, where This problem was 

solved with optimality and feasibility tolerances of 5 x 10“^ and 5 x 10“^, respectively. 
The exact solution was used as an initial guess. Figure 32 shows a comparison of the exact 
solution with the GOCM-MS (Section 5.4) numerical solutions with = 10, = 50 

and Nu = 5. 
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Figure 32: Exact solution and GOCM-MS approximation with = 10, = 50 and 

Nu = 5 for Example 7.5. 

Unlike the multi-scale Eegendre PS method outlined in Section 3.4, the GOCM-MS pro¬ 
vides a feasible solution. In fact, the maximum error of the GOCM-MS numerical solutions 
all have magnitude 0(10“"^). It is also important to point out that the required feasibility 
tolerance of 5 x 10“^ is the same order of magnitude as the largest xi spectral coefficient 
dropped from the approximation (see Eigure 13). This observation is consistent with Re¬ 
mark 3.6. Additionally, Remark 3.6 highlights the potential for a feasible solution with the 
lower bounds: 3 < 43 < and 1 < (note the associated magnitudes of the 

Eegendre spectral coefficients in Eigure 13 for xi, X2 and u). This is in fact the case: Eig¬ 
ure 33 shows a comparison of the exact solution with the GOCM-MS numerical solutions 
with Na:^ = 3, Nx 2 = 43 and = 1, solved with optimality and feasibility tolerances of 
5 X 10“^ and 5 x 10“^, respectively. 
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Figure 33: Exact solution and GOCM-MS approximation with = 3, = 43 and 

Nu = I for Example 7.5. 

The maximum errors of the new GOCM-MS approximation (with = 3, = 43 and 

Nu = 1) remain at 0(10“^) for the control, u, and 0(10“^) for the states, Xi and X 2 - Note 
that when the full-scale approach was used in Example 7.4, a GOCM-W approximation 
order of > 40 was required for the same order of accuracies. It is clear that the size of 
the NEP for this problem has been reduced significantly. The number of decision variables 
needed for the multi-scale Galerkin optimal control formulation is nearly 33% of that re¬ 
quired of the full-scale problem solved with the GOCM-W. Additionally, this multi-scale 
approach has the potential to more efficiently solve a great number of optimal control prob¬ 
lems. The larger the dimension of the multi-scale problem, the greater the potential savings 
in computational efficiency! 
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7.6. Example 7.6: Nonlinear Two-Dimensional Multi-scale Problem 
with Fixed Initial Conditions 


Consider the following nonlinear two-dimensional problem With fixed initial con¬ 
ditions given by Gong et al. [71] of minimizing the cost function 


J 


{xi — t)‘^dt, 


(7.4) 


subject to the dynamics 


Xi(t) = sin (50x1) + X2 and X2{t) = m, 


(7.5) 


and with initial conditions 


Xi(0) = 0 and ^ 2 ( 0 ) = 1. (7.6) 

The analytic solution to this problem is given by 


Xiit) = t, 

X 2 {t) = 1 — sin (50t), 
u{t) = —50cos (50t), 

obtained via Pontryagin’s maximum principle. Figures 34 and 35 show a comparison of 
the exact state and control solutions with the GOCM-W (Section 5.1) numerical solutions, 
respectively, with = 45. 
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Figure 34: Exact state solutions and GOCM-W approximation with iV = 45 for Exam 
pie 7.6. 



t 


Eigure 35: Exact control solution and GOCM-W approximation with iV = 45 for Exam 
pie 7.6. 
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It is clear from Figures 34 and 35 that problem (7.4)-(7.6) may be eonsidered multi- 
seale, where xi is the slow state and X 2 the fast state. Consider this problem now reeast in 
the form of Problem B. Here, we diseretize the slow state, xi, fast state, X 2 , and eontrol, u, 
on LGL grids, {tj}^2o and {fj}^^Q, respeetively. In order to determine 

and Nu we eonsider the speetral eoeffieients of Xi, X 2 and u given in Figure 36. 

1 

0.5 
0 

- 0.5 


0 5 10 15 20 25 30 35 40 45 50 

n 

(a) Legendre spectral coefficients for state xi. 



(b) Legendre spectral coefficients for state X2- 


0.5 

- 0 . 5 - 

5 10 15 20 25 30 35 40 45 50 

n 

(c) Legendre spectral coefficients for control u. 

Figure 36: Legendre speetral eoeffieients for xi, X 2 and u for Example 7.6. 
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If the feasibility tolerance is set to 10“^, Remark 3.6 highlights the potential for a 
feasible solution with the bounds: 1 < 39 < and 42 < Nu. This is in fact the 

case: Figures 37 and 38 show comparisons of the exact states and control with the GOCM- 
MS (Section 5.4) numerical solutions with = 1, = 39 and iV„ = 42, solved with 

optimality and feasibility tolerances of 10“'^ and 10“^, respectively. 



Figure 37: Exact state solutions and GOCM-MS approximation with = 1 and = 
39 for Example 7.6. 
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Figure 38: Exact control solution and GOCM-MS approximation with Nu = 42 for Exam¬ 
ple 7.6. 

The maximum errors of the new GOCM-MS approximation (with = 1, Nx 2 = 39 and 
Nu = 42) are 0(10“^) for the control, u, and 0(10“®) and 0(10“^) for states, xi and X 2 , 
respectively. 

7.7. Summary 

Galerkin optimal control is a versatile family of numerical formulations for solving 
optimal control problems. The examples shown in this chapter demonstrate the potential 
of this Galerkin-based family of formulations outlined in Chapters 4-6. Three particular 
highlights of Galerkin optimal control is its ability to weakly enforce problem end condi¬ 
tions, handle problems with discontinuous solutions and its potential to reduce the size of 
multi-scale problems. 

Examples 7.1, 7.2 and 7.4 demonstrate the improved accuracy provided by the 
Galerkin optimal control formulation with weak enforcement of boundary conditions over 
strong enforcement for both boundary value problems as well as problems with incomplete 
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sets of end eonditions. In partieular, the GOCM-W shows an advantage over the GOCM- 
S for low order approximations of control solutions. Examples 7.1 and 7.4 also show 
the potential advantages of the Galerkin formulations with F-LGR and LG points, respec¬ 
tively. Example 7.3 demonstrates the effectiveness of the element-based Galerkin optimal 
control formulations (such as the GOCM-DG) when employed to approximated optimal 
control problems with discontinuous controls. Lastly, Examples 7.5 and 7.6 demonstrate 
the computational efficiency in which the multi-scale Galerkin optimal control formula¬ 
tion may solve multi-scale problems, those with states and controls that evolve on different 
timescales. In contrast to the difficulties with the multi-scale Legendre PS method (see 
Section 3.4) highlighted in Chapter 3, the GOCM-MS is shown to successfully reduce the 
size of multi-scale problems. 
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CHAPTER 8: 

CONCLUSIONS AND AREAS EOR EUTURE RESEARCH 


8.1. Dissertation Summary 

This dissertation introduced and developed the theory for a Galerkin-based family 
of numerical formulations that calculate optimal trajectories by discretizing the system 
dynamics using Galerkin numerical techniques and approximate the cost function with 
Gaussian quadrature. An important result of the Galerkin formulations are that they can be 
used to prove feasibility and consistency theorems that apply to optimal control problems 
with continuous and/or piecewise continuous controls. It was shown that Galerkin optimal 
control may be formulated in a variety of ways to allow for efficiency and/or improved 
accuracy while solving a wide range of optimal control problems. 

A highlight of Galerkin optimal control is its ability to be formulated to enforce 
boundary conditions in a weak sense, imposing end conditions only up to the accuracy 
of the numerical approximation itself. The increased approximation accuracy of the weak 
boundary formulation (particularly in the approximation of control solutions) was shown on 
several linear and nonlinear problems. It was also demonstrated that the Galerkin optimal 
control formulation with weak imposition of end conditions allows for problem discretiza¬ 
tions with other than LGL points. Galerkin optimal control with Legendre-Gauss-Radau 
and Legendre-Gauss points were shown to be advantageous due to the increased accuracy 
of solutions. Galerkin optimal control may also be formulated with other than Lagrangian 
test functions, such as the Legendre polynomials. 

In addition, Galerkin optimal control has proven to be effective in reducing the 
dimension of multi-scale problems, those in which states and controls evolve on different 
timescales. In one example presented the number of decision variables required by the 
multi-scale Galerkin optimal control formulation was nearly 33% of that required of the 
full-scale problem. The multi-scale formulation has the potential to more efficiently solve 
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a great number of optimal control problems, certainly those that include fast and slow 
dynamics. 

For optimal control problems with discontinuities (such as bang-bang control prob¬ 
lems), an element-based approach has shown to be beneficial. In general, using the element- 
based Galerkin optimal control formations discussed (both continuous and discontinuous 
Galerkin) may lead to higher computational efficiencies and the formulation may be re¬ 
tooled to incorporate hp-adaptive techniques, such as the spectral algorithm discussed in 
[86]. Additionally, the discontinuous Galerkin formulation may be advantageous from a 
parallel computing standpoint. 

Galerkin optimal control has demonstrated exponential convergence for a large 
class of problems. It is clear that Galerkin optimal control is a versatile and accurate family 
of formulations that has the potential to provide real time optimal control solutions for a 
number of applications. 

8.2. Future Work 

Galerkin optimal control shows the potential for solving a wide range of optimal 
control problems with a variety of state and control constraints. However, application of 
the Galerkin optimal control formulations to real-world problems have been somewhat 
limited due to research time limitations. Future application of Galerkin optimal control to 
problems with real-world conditions is necessary to demonstrate its true versitiliy. Testing 
Galerkin optimal control on different types of control problems will inevitably highlight the 
strengths (and weaknesses) of each formulation. Lastly, this dissertation includes a number 
of important theorems that serve as the theoretical foundations for Galerkin optimal control, 
however, the list is not complete. Future work to increase the theoretical underpinnings of 
Galerkin optimal control, to include a rate of convergence analysis, is forthcoming. 
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APPENDIX A: 

NORMS AND FUNCTIONAL SPACES 


Throughout this dissertation, the spaces, for 1 < p < cx), are used quite often. 
They consist of all measurable functions n : [—1,1] —)■ M, such that [87] 


\v{t)fdt < oo, 


'-1 


with finite L^-norm, 






\v{t)\^dt < oo 


In particular, the space is referenced quite readily, defined by 


>1 


\v\\^2 = ( I Ht)\ dt ) , 


which is induced by the inner product 

{u,v) = j u{t)v{t)dt. 

Additionally, the space consist of all measurable functions v : [—1,1] —)■ 
such that |n(t) I < M for almost alH G [—1,1]. The L°°-norm can be expressed as 

ll'yll^oo = inf{M||n(t)| < M almost everywhere on t G [—1,1]}. 


The idea of a function of bounded variation is also used within this dissertation. To 
frame this idea first we must define the total variation, V(u). For a function u : [—1,1] —)■ M 
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the total variation of m on [—1,1] is defined as 

N 

V{u) = sup{^|M(tj) - u{tj.i)\ {tk}k=o e P}, 
i=i 

where P is the set of all finite partitions of [-1,1]. If the total variation is bounded, V(u) < 
oo, then u is said to be of bounded variation in [—1,1], 
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APPENDIX B: 
DISCRETE NORMS 


Let grid be associated with quadrature weights Section 2.2.2). 

Then the discrete norm, ||n||^, from [46] is defined as the quantity 

1 


where the discrete inner product is accomplished by numerical integration and given by 


N 

{v,v)n = 

j=0 


In the case that v'^ G P 2 N+S and LG, LGR or LGL points, the numerical 

integration is exact and 


{v,v)n 



u^dt, 


where S = 1,0, —1 for LG, LGR or LGL points, respectively, and are the associ¬ 

ated weights. 

Additionally, there exist constants a,/3 > 0 such that 


Oi\\v\\^2 < ||n||^ < /3\\v\\^2 


for all V e Pn- 

Lastly, ||n||g^ denotes the maximum element of vector, v eMP. 
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APPENDIX C: 
SOBOLEV SPACES 


Throughout this dissertation, the Sobolev spaces, are referenced quite often. 

They consist of all functions, v : [—1,1] —)■ having weak derivative, G L^, where 

0 < i < m, with the norm [87] 


I 



J=0 





where ||n||ip denotes. 


LP ~ 


\v(t)\^dt 


'-1 


1 

V 


The seminorm may be expressed as 


V 


y^m,p-,N 


1 



Additionally, Sobolev spaces may be defined with a fractional order. They consist of all 
measurable functions v : [—1,1] —)■ M, such that [88] 


W--<- = {p e L' : / ' < oo}, 


'-1 J-I 


\x - y\ 


for 0 < a < 1 and 1 < p < cx), with the norm. 


\W^''P 


p \v{x) — v{y)f 

\v{t)\^dt+ I I - - ^T^J-dxdy 




\x - y\ 


Lastly, due to the extensive use of the space notation is simplified by letting = 


H^. 
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